In [1]:
rm(list=ls())
gc()
setwd("/hpc/group/pbenfeylab/CheWei/")
A matrix: 2 x 6 of type dbl
used(Mb)gc trigger(Mb)max used(Mb)
Ncells 62527533.4136147072.8120957264.6
Vcells1159072 8.9838860864.0180227913.8
In [2]:
as.numeric(system("awk '/MemFree/ {print $2}' /proc/meminfo", intern=TRUE))
196300512
In [3]:
suppressMessages(library(Seurat))
suppressMessages(library(tricycle))
suppressMessages(library(cowplot))
suppressMessages(library(ggplot2))
suppressMessages(library(scattermore))
suppressMessages(library(scater))
suppressMessages(library(cowplot))
suppressMessages(library(RColorBrewer))
suppressMessages(library(grid))
suppressMessages(library(gplots))
suppressMessages(library(tidyverse))

prep <- function(seu){
seu$time.anno.Li.crude <- gsub("_.*$","",seu$time.celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("Distal Columella","Distal Columella_Columella",seu$time.celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("Distal Lateral Root Cap","Distal Lateral Root Cap_Lateral Root Cap",seu$celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("Proximal Columella","Proximal Columella_Columella",seu$celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("Proximal Lateral Root Cap","Proximal Lateral Root Cap_Lateral Root Cap",seu$celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("^.*_","",seu$celltype.anno.Li.crude)
    return(seu)
}
plot_anno <- function(rc.integrated){
order <- c("Quiescent Center", "Ground Tissue","Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem","Protophloem", "Xylem", "Procambium","Pericycle","Phloem Pole Pericycle", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#dd77ec", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$celltype.anno.Li.crude <- factor(rc.integrated$celltype.anno.Li.crude, levels = order[sort(match(unique(rc.integrated$celltype.anno.Li.crude),order))]) 
color <- palette[sort(match(unique(rc.integrated$celltype.anno.Li.crude),order))]
p1 <- DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno.Li.crude", cols=color)
p2 <- DimPlot(rc.integrated, reduction = "umap", group.by = "time.anno.Li.crude", order = c("Distal Columella","Proximal Columella","Distal Lateral Root Cap","Proximal Lateral Root Cap","Maturation","Elongation", "Transition Domain", "Proliferation Domain"),
        cols = c('#F7E7B0','#FFC400','#2B871F','#005E3B',"#deebf7", "#3182bd", '#fee0d2','#de2d26'))
p3 <- DimPlot(rc.integrated, reduction = "umap", group.by = "consensus.time.group", cols=brewer.pal(10,"Spectral"))
p4 <- DimPlot(rc.integrated, reduction = "umap", group.by = "branch.anno")
options(repr.plot.width=16, repr.plot.height=12)
gl <- lapply(list(p1, p2, p3, p4), ggplotGrob)
gwidth <- do.call(unit.pmax, lapply(gl, "[[", "widths"))
gl <- lapply(gl, "[[<-", "widths", value = gwidth)
gridExtra::grid.arrange(grobs=gl, ncol=2)
}
In [4]:
sessionInfo()
R version 4.2.2 (2022-10-31)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: AlmaLinux 9.3 (Shamrock Pampas Cat)

Matrix products: default
BLAS/LAPACK: /hpc/group/pbenfeylab/ch416/miniconda3/envs/seu4/lib/libopenblasp-r0.3.21.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] grid      stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] forcats_0.5.2               stringr_1.5.1              
 [3] dplyr_1.1.3                 purrr_1.0.2                
 [5] readr_2.1.3                 tidyr_1.3.0                
 [7] tibble_3.2.1                tidyverse_1.3.2            
 [9] gplots_3.1.3                RColorBrewer_1.1-3         
[11] scater_1.26.1               scuttle_1.8.0              
[13] scattermore_1.2             ggplot2_3.4.4              
[15] cowplot_1.1.1               tricycle_1.6.0             
[17] SingleCellExperiment_1.20.0 SummarizedExperiment_1.28.0
[19] Biobase_2.58.0              GenomicRanges_1.50.0       
[21] GenomeInfoDb_1.34.8         IRanges_2.32.0             
[23] S4Vectors_0.36.0            BiocGenerics_0.44.0        
[25] MatrixGenerics_1.10.0       matrixStats_1.1.0          
[27] SeuratObject_4.1.3          Seurat_4.1.1.9001          

loaded via a namespace (and not attached):
  [1] utf8_1.2.4                spatstat.explore_3.2-5   
  [3] reticulate_1.34.0         tidyselect_1.2.0         
  [5] RSQLite_2.3.1             AnnotationDbi_1.60.0     
  [7] htmlwidgets_1.6.2         BiocParallel_1.32.5      
  [9] Rtsne_0.16                munsell_0.5.0            
 [11] ScaledMatrix_1.6.0        codetools_0.2-19         
 [13] ica_1.0-3                 pbdZMQ_0.3-8             
 [15] future_1.33.0             miniUI_0.1.1.1           
 [17] withr_2.5.2               spatstat.random_3.2-1    
 [19] colorspace_2.1-0          progressr_0.14.0         
 [21] uuid_1.1-0                ROCR_1.0-11              
 [23] tensor_1.5                listenv_0.9.0            
 [25] repr_1.1.4                GenomeInfoDbData_1.2.9   
 [27] polyclip_1.10-6           bit64_4.0.5              
 [29] parallelly_1.36.0         vctrs_0.6.4              
 [31] generics_0.1.3            circular_0.4-95          
 [33] timechange_0.1.1          R6_2.5.1                 
 [35] ggbeeswarm_0.7.1          rsvd_1.0.5               
 [37] bitops_1.0-7              spatstat.utils_3.0-4     
 [39] cachem_1.0.8              DelayedArray_0.24.0      
 [41] assertthat_0.2.1          promises_1.2.1           
 [43] scales_1.2.1              googlesheets4_1.0.1      
 [45] beeswarm_0.4.0            gtable_0.3.4             
 [47] beachmat_2.14.0           globals_0.16.2           
 [49] goftest_1.2-3             rlang_1.1.2              
 [51] splines_4.2.2             lazyeval_0.2.2           
 [53] gargle_1.2.1              broom_1.0.2              
 [55] spatstat.geom_3.2-7       modelr_0.1.10            
 [57] reshape2_1.4.4            abind_1.4-5              
 [59] backports_1.4.1           httpuv_1.6.12            
 [61] tools_4.2.2               ellipsis_0.3.2           
 [63] ggridges_0.5.4            Rcpp_1.0.11              
 [65] plyr_1.8.9                base64enc_0.1-3          
 [67] sparseMatrixStats_1.10.0  zlibbioc_1.44.0          
 [69] RCurl_1.98-1.6            deldir_1.0-9             
 [71] pbapply_1.7-2             viridis_0.6.4            
 [73] zoo_1.8-12                haven_2.5.1              
 [75] ggrepel_0.9.4             cluster_2.1.4            
 [77] fs_1.6.3                  magrittr_2.0.3           
 [79] data.table_1.14.8         RSpectra_0.16-1          
 [81] reprex_2.0.2              lmtest_0.9-40            
 [83] RANN_2.6.1                googledrive_2.0.0        
 [85] mvtnorm_1.1-3             ggnewscale_0.4.8         
 [87] fitdistrplus_1.1-11       hms_1.1.2                
 [89] patchwork_1.1.3           mime_0.12                
 [91] evaluate_0.23             xtable_1.8-4             
 [93] readxl_1.4.1              fastDummies_1.7.3        
 [95] gridExtra_2.3             compiler_4.2.2           
 [97] KernSmooth_2.23-20        crayon_1.5.2             
 [99] htmltools_0.5.7           tzdb_0.3.0               
[101] later_1.3.1               lubridate_1.9.0          
[103] DBI_1.1.3                 dbplyr_2.2.1             
[105] MASS_7.3-58.3             boot_1.3-28.1            
[107] Matrix_1.5-4              cli_3.6.1                
[109] parallel_4.2.2            igraph_1.5.1             
[111] pkgconfig_2.0.3           sp_2.1-1                 
[113] IRdisplay_1.1             plotly_4.10.3            
[115] spatstat.sparse_3.0-3     xml2_1.3.3               
[117] vipor_0.4.5               XVector_0.38.0           
[119] rvest_1.0.3               digest_0.6.33            
[121] sctransform_0.4.1         RcppAnnoy_0.0.21         
[123] spatstat.data_3.0-3       Biostrings_2.66.0        
[125] cellranger_1.1.0          leiden_0.4.3             
[127] uwot_0.1.16               DelayedMatrixStats_1.20.0
[129] shiny_1.7.5.1             gtools_3.9.4             
[131] lifecycle_1.0.4           nlme_3.1-162             
[133] jsonlite_1.8.7            BiocNeighbors_1.16.0     
[135] viridisLite_0.4.2         fansi_1.0.5              
[137] pillar_1.9.0              lattice_0.21-8           
[139] KEGGREST_1.38.0           fastmap_1.1.1            
[141] httr_1.4.7                survival_3.4-0           
[143] glue_1.6.2                png_0.1-8                
[145] bit_4.0.5                 stringi_1.8.1            
[147] blob_1.2.4                RcppHNSW_0.5.0           
[149] BiocSingular_1.14.0       caTools_1.18.2           
[151] memoise_2.0.1             IRkernel_1.3.1.9000      
[153] irlba_2.3.5.1             future.apply_1.11.0      
In [5]:
cc.genes.goldy <- read.csv("./tricycle/Goldy_2021_CycB1-1-sorted_Enriched_Depleted_G2M_rm-proto-genes.csv")
cc.genes.zhang <- read.csv("./tricycle/Zhang_2021_Dev_cell_Core_cell_cycle_genes.csv")
cc.genes.amigo.go <- as.character(read.table("./tricycle/GO_cell_cycle_annotation_arabidopsis.txt", header=F)$V1)
cc.genes.f <- unique(c(cc.genes.zhang$Gene, cc.genes.goldy$gene, cc.genes.amigo.go))
length(cc.genes.f)
pp.genes <- as.character(read.table("./CW_data/escoring/Root_sc/Protoplasting_DEgene_FC2_list.txt", header=F)$V1)
3241
In [6]:
table(cc.genes.zhang$Phase)
G1/G0    G2     M     S 
   18    95    77    53 
In [7]:
head(cc.genes.zhang)
A data.frame: 6 x 5
GenePhaseSymbolAraport11Function
<chr><chr><chr><chr><chr>
1AT1G12845G1/G0NA transmembrane protein;(source:Araport11) NA
2AT3G26960G1/G0NA Pollen Ole e 1 allergen and extensin family protein;(source:Araport11) NA
3AT5G02260G1/G0ATEXP9 expansin A9;(source:Araport11) member of Alpha-Expansin Gene Family. Naming convention from the Expansin Working Group (Kende et al, 2004. Plant Mol Bio)
4AT1G33930G1/G0NA P-loop containing nucleoside triphosphate hydrolases superfamily protein;(source:Araport11)NA
5AT5G01445G1/G0NA hypothetical protein;(source:Araport11) NA
6AT5G67260G1/G0CYCD3;2CYCLIN D3;(source:Araport11) Encode CYCD3;2, a CYCD3 D-type cyclin. Important for determining cell number in developing lateral organs. Mediating cytokinin effects in apical growth and development.
In [8]:
## Cell Cyle genes used
cc.gene.list <- as.list(cc.genes.f)

Load reference¶

In [9]:
rc.integrated <- readRDS("./scRNA-seq/Integrated_Objects/rc.integrated_18S_WT_cell_cycle_seu4_245genes_Goldy_Zhang_AMIGO_20230825.rds")
In [10]:
rownames(rc.integrated)
  1. 'AT1G02730'
  2. 'AT5G56120'
  3. 'AT5G13840'
  4. 'AT1G63100'
  5. 'AT1G03780'
  6. 'AT5G15510'
  7. 'AT1G16630'
  8. 'AT3G51230'
  9. 'AT4G35620'
  10. 'AT4G02800'
  11. 'AT4G05190'
  12. 'AT3G55660'
  13. 'AT3G11520'
  14. 'AT1G20610'
  15. 'AT1G18370'
  16. 'AT1G76310'
  17. 'AT4G14330'
  18. 'AT5G48310'
  19. 'AT5G55520'
  20. 'AT4G01730'
  21. 'AT1G44110'
  22. 'AT5G45700'
  23. 'AT4G38062'
  24. 'AT3G56100'
  25. 'AT4G26660'
  26. 'AT5G60930'
  27. 'AT4G09060'
  28. 'AT4G37490'
  29. 'AT5G51600'
  30. 'AT3G23670'
  31. 'AT3G44050'
  32. 'AT1G11220'
  33. 'AT2G22610'
  34. 'AT1G04425'
  35. 'AT4G11080'
  36. 'AT4G17000'
  37. 'AT3G19050'
  38. 'AT1G53140'
  39. 'AT3G15550'
  40. 'AT1G72670'
  41. 'AT1G34355'
  42. 'AT1G05440'
  43. 'AT2G44190'
  44. 'AT3G54710'
  45. 'AT5G67270'
  46. 'AT4G21270'
  47. 'AT1G72250'
  48. 'AT3G58650'
  49. 'AT1G61450'
  50. 'AT1G23790'
  51. 'AT5G55830'
  52. 'AT5G44040'
  53. 'AT1G59540'
  54. 'AT5G09500'
  55. 'AT5G43080'
  56. 'AT4G18570'
  57. 'AT5G23910'
  58. 'AT3G57540'
  59. 'AT1G69770'
  60. 'AT4G28430'
  61. 'AT5G47500'
  62. 'AT1G75640'
  63. 'AT1G76740'
  64. 'AT4G30860'
  65. 'AT3G06840'
  66. 'AT2G36200'
  67. 'AT3G25100'
  68. 'AT1G33670'
  69. 'AT4G05520'
  70. 'AT5G60150'
  71. 'AT4G21070'
  72. 'AT3G26050'
  73. 'AT5G52220'
  74. 'AT2G44580'
  75. 'AT2G47920'
  76. 'AT5G46740'
  77. 'AT3G27640'
  78. 'AT4G14150'
  79. 'AT3G12170'
  80. 'AT5G42700'
  81. 'AT3G48490'
  82. 'AT3G24495'
  83. 'AT1G06310'
  84. 'AT5G01420'
  85. 'AT1G10780'
  86. 'AT3G52115'
  87. 'AT2G37420'
  88. 'AT1G10400'
  89. 'AT1G80190'
  90. 'AT4G00020'
  91. 'AT2G38620'
  92. 'AT1G77630'
  93. 'AT4G21820'
  94. 'AT1G29418'
  95. 'AT5G55820'
  96. 'AT2G42260'
  97. 'AT3G23740'
  98. 'AT5G11780'
  99. 'AT1G26330'
  100. 'AT2G32590'
  101. 'AT3G54750'
  102. 'AT4G39590'
  103. 'AT4G17240'
  104. 'AT1G50240'
  105. 'AT5G59240'
  106. 'AT3G59550'
  107. 'AT5G11510'
  108. 'AT3G05740'
  109. 'AT1G07270'
  110. 'AT3G55340'
  111. 'AT1G20720'
  112. 'AT1G75150'
  113. 'AT3G61150'
  114. 'AT3G58100'
  115. 'AT2G20980'
  116. 'AT4G39630'
  117. 'AT3G18680'
  118. 'AT5G62170'
  119. 'AT3G09080'
  120. 'AT1G63650'
  121. 'AT2G02680'
  122. 'AT4G22970'
  123. 'AT5G20850'
  124. 'AT3G48160'
  125. 'AT3G18730'
  126. 'AT5G14660'
  127. 'AT5G25380'
  128. 'AT2G29680'
  129. 'AT2G42290'
  130. 'AT5G56220'
  131. 'AT3G63480'
  132. 'AT5G07810'
  133. 'AT5G18570'
  134. 'AT5G62410'
  135. 'AT1G77550'
  136. 'AT2G24970'
  137. 'AT5G40840'
  138. 'AT4G13370'
  139. 'AT1G09700'
  140. 'AT1G72480'
  141. 'AT1G53542'
  142. 'AT4G02150'
  143. 'AT3G61610'
  144. 'AT3G59430'
  145. 'AT2G04530'
  146. 'AT5G52170'
  147. 'AT3G62110'
  148. 'AT5G05900'
  149. 'AT2G40570'
  150. 'AT4G34960'
  151. 'AT1G19520'
  152. 'AT2G34450'
  153. 'AT2G45590'
  154. 'AT1G30960'
  155. 'AT1G54960'
  156. 'AT1G72820'
  157. 'AT1G05520'
  158. 'AT5G37630'
  159. 'AT1G11125'
  160. 'AT2G37560'
  161. 'AT3G03810'
  162. 'AT3G24430'
  163. 'AT3G19720'
  164. 'AT4G23510'
  165. 'AT1G12800'
  166. 'AT4G12540'
  167. 'AT5G13830'
  168. 'AT4G12740'
  169. 'AT4G16700'
  170. 'AT5G47400'
  171. 'AT3G13150'
  172. 'AT5G63920'
  173. 'AT5G66770'
  174. 'AT3G09410'
  175. 'AT1G71060'
  176. 'AT4G38230'
  177. 'AT5G01470'
  178. 'AT2G42190'
  179. 'AT5G51560'
  180. 'AT4G14980'
  181. 'AT4G29170'
  182. 'AT1G20410'
  183. 'AT4G31400'
  184. 'AT3G46020'
  185. 'AT1G10300'
  186. 'AT2G17670'
  187. 'AT5G35520'
  188. 'AT1G71760'
  189. 'AT1G73180'
  190. 'AT3G13060'
  191. 'AT5G59740'
  192. 'AT4G34360'
  193. 'AT2G20810'
  194. 'AT1G49910'
  195. 'AT1G04710'
  196. 'AT1G75140'
  197. 'AT1G04130'
  198. 'AT4G12620'
  199. 'AT2G21800'
  200. 'AT4G01810'
  201. 'AT2G04660'
  202. 'AT5G60840'
  203. 'AT1G79740'
  204. 'AT5G38880'
  205. 'AT1G61000'
  206. 'AT4G15890'
  207. 'AT3G62620'
  208. 'AT2G46380'
  209. 'AT5G12440'
  210. 'AT4G19710'
  211. 'AT5G06590'
  212. 'AT3G60740'
  213. 'AT5G22340'
  214. 'AT5G63120'
  215. 'AT4G24175'
  216. 'AT3G07630'
  217. 'AT4G28020'
  218. 'AT1G21150'
  219. 'AT1G47230'
  220. 'AT5G02820'
  221. 'AT4G17430'
  222. 'AT4G19130'
  223. 'AT5G63690'
  224. 'AT3G20780'
  225. 'AT1G27060'
  226. 'AT2G04940'
  227. 'AT1G62150'
  228. 'AT3G19570'
  229. 'AT3G58580'
  230. 'AT5G09380'
  231. 'AT1G09280'
  232. 'AT1G67780'
  233. 'AT2G33610'
  234. 'AT4G19185'
  235. 'AT5G46850'
  236. 'AT5G12360'
  237. 'AT2G03670'
  238. 'AT5G03415'
  239. 'AT5G06180'
  240. 'AT4G04190'
  241. 'AT5G45760'
  242. 'AT2G37680'
  243. 'AT4G15850'
  244. 'AT4G26720'
  245. 'AT1G07020'
In [11]:
length(rownames(rc.integrated))
245
In [12]:
write.csv(data.frame("GeneID"=rownames(rc.integrated)),"./tradeseq/245_cell_cycle_related_genes_for_reference.csv")
In [10]:
## Prepare New Ref (All cells)
corrected.m <- rc.integrated@assays$integrated@data
In [12]:
# run PCA on the batch effects corrected matrix and get the rotaions scores for the top 2 PCs
pca.m <- scater::calculatePCA(corrected.m, ntop = 500)
In [13]:
new.ref <- attr(pca.m, 'rotation')[, seq_len(2)]
In [15]:
saveRDS(new.ref, "./tradeseq/Tricycle_Reference_Arabidopsis_Root.rds")
In [14]:
unique(rc.integrated$orig.ident)
  1. 'dc1'
  2. 'dc2'
  3. 'sc_10_at'
  4. 'sc_11'
  5. 'sc_12'
  6. 'tnw1'
  7. 'sc_20'
  8. 'sc_71'
  9. 'sc_157'
  10. 'sc_158'
  11. 'sc_159'
  12. 'sc_169'
  13. 'sc_170'
  14. 'sc_171'
  15. 'sc_215'
  16. 'sc_221'
  17. 'sc_229'
  18. 'sc_235'
In [13]:
#suppressMessages(rc.integrated <- RunUMAP(rc.integrated, reduction = "pca", dims = 1:8))
In [14]:
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "umap", group.by = "seurat_clusters", label = TRUE)+theme(legend.position="none")
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [15]:
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "pca", group.by = "seurat_clusters", label = TRUE)+theme(legend.position="none")
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [16]:
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "umap", group.by = "orig.ident", label = FALSE, pt.size=1, order=FALSE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [17]:
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "pca", group.by = "orig.ident", label = FALSE, pt.size=1, order=FALSE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [18]:
order <- c("Quiescent Center", "Ground Tissue","Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem","Protophloem", "Xylem", "Procambium","Xylem Pole Pericycle","Phloem Pole Pericycle", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#dd77ec", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$celltype.anno.Li <- factor(rc.integrated$celltype.anno.Li, levels = order[sort(match(unique(rc.integrated$celltype.anno.Li),order))]) 
color <- palette[sort(match(unique(rc.integrated$celltype.anno.Li),order))]
In [19]:
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno.Li", label = FALSE, cols=color, pt.size=1, order=FALSE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [20]:
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "pca", group.by = "celltype.anno.Li", label = FALSE, cols=color, pt.size=1, order=FALSE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [21]:
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "pca", group.by = "time.anno.Li", label = FALSE, pt.size=1, order=FALSE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [22]:
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "pca", group.by = "consensus.time.group", label = FALSE, pt.size=1, order=FALSE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [23]:
## S-phase cluster
sort(table(rc.integrated[,which(rc.integrated$seurat_clusters==0)]$celltype.anno.Li))
     Quiescent Center                Phloem Phloem Pole Pericycle 
                   38                   275                   323 
            Columella                 Xylem                Cortex 
                  331                   795                  1106 
           Procambium  Xylem Pole Pericycle            Endodermis 
                 1426                  1439                  1592 
         Atrichoblast           Trichoblast      Lateral Root Cap 
                 1646                  2546                  2879 
In [24]:
## S-phase cluster
sort(table(rc.integrated[,which(rc.integrated$seurat_clusters==0)]$time.anno.Li))
  Distal Lateral Root Cap          Distal Columella        Proximal Columella 
                      136                       161                       170 
               Maturation         Transition Domain Proximal Lateral Root Cap 
                      763                      1900                      2743 
     Proliferation Domain                Elongation 
                     3610                      4913 
In [25]:
## G2M-phase cluster
rev(sort(table(rc.integrated[,which(rc.integrated$seurat_clusters==4)]$celltype.anno.Li)))
     Lateral Root Cap          Atrichoblast            Endodermis 
                 2981                  1013                   790 
          Trichoblast                Cortex                 Xylem 
                  550                   429                   284 
            Columella  Xylem Pole Pericycle            Procambium 
                  275                   252                   249 
               Phloem Phloem Pole Pericycle      Quiescent Center 
                  233                   100                     8 
In [26]:
## G2M-phase cluster
rev(sort(table(rc.integrated[,which(rc.integrated$seurat_clusters==4)]$time.anno.Li)))
Proximal Lateral Root Cap      Proliferation Domain                Elongation 
                     2863                      1840                      1324 
               Maturation         Transition Domain          Distal Columella 
                      533                       211                       185 
  Distal Lateral Root Cap        Proximal Columella 
                      118                        90 
In [27]:
## G2M-phase cluster
rev(sort(table(rc.integrated[,which(rc.integrated$seurat_clusters==4)]$time.celltype.anno.Li)))
                       Proximal Lateral Root Cap 
                                            2863 
               Proliferation Domain_Atrichoblast 
                                             828 
                 Proliferation Domain_Endodermis 
                                             509 
                     Proliferation Domain_Cortex 
                                             264 
                Proliferation Domain_Trichoblast 
                                             210 
                 Elongation_Xylem Pole Pericycle 
                                             196 
                           Elongation_Endodermis 
                                             195 
                                Distal Columella 
                                             185 
                           Elongation_Procambium 
                                             176 
          Elongation_Metaphloem & Companion Cell 
                                             153 
                   Transition Domain_Trichoblast 
                                             144 
                           Elongation_Protoxylem 
                                             128 
                           Maturation_Protoxylem 
                                             123 
                               Elongation_Cortex 
                                             118 
                         Distal Lateral Root Cap 
                                             118 
                          Elongation_Trichoblast 
                                             114 
                         Elongation_Atrichoblast 
                                              99 
                              Proximal Columella 
                                              90 
                Elongation_Phloem Pole Pericycle 
                                              90 
                           Maturation_Endodermis 
                                              85 
                          Maturation_Trichoblast 
                                              82 
                           Maturation_Procambium 
                                              71 
                         Maturation_Atrichoblast 
                                              57 
                 Maturation_Xylem Pole Pericycle 
                                              41 
                               Maturation_Cortex 
                                              33 
                  Transition Domain_Atrichoblast 
                                              29 
                          Elongation_Protophloem 
                                              29 
          Maturation_Metaphloem & Companion Cell 
                                              28 
                            Elongation_Metaxylem 
                                              26 
   Transition Domain_Metaphloem & Companion Cell 
                                              14 
                        Transition Domain_Cortex 
                                              14 
       Proliferation Domain_Xylem Pole Pericycle 
                                              12 
           Proliferation Domain_Quiescent Center 
                                               8 
                Maturation_Phloem Pole Pericycle 
                                               8 
                  Proliferation Domain_Metaxylem 
                                               5 
                          Maturation_Protophloem 
                                               5 
          Transition Domain_Xylem Pole Pericycle 
                                               3 
                   Transition Domain_Protophloem 
                                               2 
         Transition Domain_Phloem Pole Pericycle 
                                               2 
Proliferation Domain_Metaphloem & Companion Cell 
                                               2 
                    Transition Domain_Procambium 
                                               1 
                     Transition Domain_Metaxylem 
                                               1 
                    Transition Domain_Endodermis 
                                               1 
                 Proliferation Domain_Protoxylem 
                                               1 
                 Proliferation Domain_Procambium 
                                               1 

Top Genes contribute the most in each PCs¶

In [30]:
pc1r <- as.data.frame(rc.integrated[['pca']][,1]) %>% arrange(desc(PC_1)) %>% rownames(.)
In [31]:
pc1r[1:9]
  1. 'AT1G18370'
  2. 'AT1G02730'
  3. 'AT5G15510'
  4. 'AT1G16630'
  5. 'AT5G56120'
  6. 'AT3G55660'
  7. 'AT4G14330'
  8. 'AT4G05190'
  9. 'AT3G19050'
In [32]:
## PC1 positive
FeaturePlot(rc.integrated, reduction="pca", features = c(pc1r[1:9]), cols = c("#EEEEEE","#FF4040"), order=TRUE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [33]:
## PC1 positive
FeaturePlot(rc.integrated, reduction="umap", features = c(pc1r[1:9]), cols = c("#EEEEEE","#FF4040"), order=TRUE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [34]:
pc1r[1:20]
  1. 'AT1G18370'
  2. 'AT1G02730'
  3. 'AT5G15510'
  4. 'AT1G16630'
  5. 'AT5G56120'
  6. 'AT3G55660'
  7. 'AT4G14330'
  8. 'AT4G05190'
  9. 'AT3G19050'
  10. 'AT4G26660'
  11. 'AT1G63100'
  12. 'AT1G03780'
  13. 'AT5G13840'
  14. 'AT5G60930'
  15. 'AT5G55520'
  16. 'AT3G44050'
  17. 'AT1G20610'
  18. 'AT4G02800'
  19. 'AT4G09060'
  20. 'AT4G38062'
In [35]:
## PC1 negative
FeaturePlot(rc.integrated, reduction="pca", features = c(rev(pc1r)[1:9]), cols = c("#EEEEEE","#FF4040"),order=TRUE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [36]:
pc2r <- as.data.frame(rc.integrated[['pca']][,2]) %>% arrange(desc(PC_2)) %>% rownames(.)
In [37]:
pc2r[1:20]
  1. 'AT5G56120'
  2. 'AT5G15510'
  3. 'AT1G02730'
  4. 'AT4G26660'
  5. 'AT1G03780'
  6. 'AT3G55660'
  7. 'AT1G16630'
  8. 'AT4G14330'
  9. 'AT5G55520'
  10. 'AT5G13840'
  11. 'AT4G02800'
  12. 'AT2G44190'
  13. 'AT5G44040'
  14. 'AT1G63100'
  15. 'AT3G51230'
  16. 'AT3G11520'
  17. 'AT4G28430'
  18. 'AT3G15550'
  19. 'AT4G17000'
  20. 'AT4G09060'
In [38]:
rev(pc2r)[1:20]
  1. 'AT1G69770'
  2. 'AT3G25100'
  3. 'AT5G43080'
  4. 'AT4G00020'
  5. 'AT3G05740'
  6. 'AT3G54750'
  7. 'AT3G27640'
  8. 'AT4G30860'
  9. 'AT1G80190'
  10. 'AT3G52115'
  11. 'AT5G46740'
  12. 'AT4G21070'
  13. 'AT2G44580'
  14. 'AT3G12170'
  15. 'AT2G29680'
  16. 'AT5G20850'
  17. 'AT3G24495'
  18. 'AT5G52220'
  19. 'AT3G18730'
  20. 'AT2G42260'
In [39]:
## PC2 positive
FeaturePlot(rc.integrated, reduction="pca", features = c(pc2r[1:9]), cols = c("#EEEEEE","#FF4040"),order=TRUE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [40]:
## PC2 negative
FeaturePlot(rc.integrated, reduction="pca", features = c(rev(pc2r)[1:9]), cols = c("#EEEEEE","#FF4040"),order=TRUE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [41]:
## PC2 negative
FeaturePlot(rc.integrated, reduction="umap", features = c(rev(pc2r)[1:9]), cols = c("#EEEEEE","#FF4040"),order=TRUE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [42]:
pc3r <- as.data.frame(rc.integrated[['pca']][,3]) %>% arrange(desc(PC_3)) %>% rownames(.)
In [43]:
pc3r[1:20]
  1. 'AT1G29418'
  2. 'AT5G59240'
  3. 'AT3G55340'
  4. 'AT5G60930'
  5. 'AT1G73180'
  6. 'AT1G77550'
  7. 'AT4G18570'
  8. 'AT1G19520'
  9. 'AT3G07630'
  10. 'AT1G72670'
  11. 'AT3G13150'
  12. 'AT1G30960'
  13. 'AT1G53542'
  14. 'AT5G62410'
  15. 'AT1G09700'
  16. 'AT1G53140'
  17. 'AT1G04710'
  18. 'AT2G36200'
  19. 'AT3G19050'
  20. 'AT1G04130'
In [44]:
rev(pc3r)[1:20]
  1. 'AT3G25100'
  2. 'AT5G43080'
  3. 'AT5G56120'
  4. 'AT3G05740'
  5. 'AT3G52115'
  6. 'AT4G00020'
  7. 'AT3G27640'
  8. 'AT3G51230'
  9. 'AT4G30860'
  10. 'AT2G29680'
  11. 'AT1G03780'
  12. 'AT1G63100'
  13. 'AT1G02730'
  14. 'AT4G02800'
  15. 'AT5G15510'
  16. 'AT3G18730'
  17. 'AT4G26660'
  18. 'AT4G21070'
  19. 'AT4G17000'
  20. 'AT5G46740'

Phase signature¶

In [45]:
zscore <- function(x){(x-mean(x))/sd(x)}

Zhang's list¶

In [46]:
cc.genes.zhang <- cc.genes.zhang[match(intersect(rownames(rc.integrated),cc.genes.zhang$Gene),cc.genes.zhang$Gene),] %>% arrange(Phase)
In [47]:
head(cc.genes.zhang)
A data.frame: 6 x 5
GenePhaseSymbolAraport11Function
<chr><chr><chr><chr><chr>
1AT1G16630G2NA transmembrane protein;(source:Araport11) NA
2AT4G05190G2ATK5 kinesin 5;(source:Araport11) ATK5 encodes a kinesin protein involved in microtubule spindle morphogenesis. It acts as a minus-end directed motor as well as a plus-end tracking protein (+TIP). Localizes to mitotic spindle midzones and regions rich in growing plus-ends within phragmoplasts.
3AT1G20610G2CYCB2;3Cyclin B2;(source:Araport11) NA
4AT5G55520G2NA kinesin-like protein;(source:Araport11) NA
5AT4G01730G2NA DHHC-type zinc finger family protein;(source:Araport11)NA
6AT1G44110G2CYCA1;1Cyclin A1;(source:Araport11) NA
In [17]:
head(cc.genes.zhang %>% filter(Phase == "S"))
A data.frame: 6 x 5
GenePhaseSymbolAraport11Function
<chr><chr><chr><chr><chr>
1AT5G10400SH3.1 Histone superfamily protein;(source:Araport11)NA
2AT5G22880SH2B histone B2;(source:Araport11) Encodes a histone 2B (H2B) protein. This protein can be ubiquitinated in planta, and this modification depends on the HUB1 and HUB2 E3 ubiquitin ligases.
3AT5G59870Sh2a.w.6histone H2A 6;(source:Araport11) Encodes HTA6, a histone H2A protein.
4AT5G59690SNA Histone superfamily protein;(source:Araport11)NA
5AT5G65360SH3.1 Histone superfamily protein;(source:Araport11)NA
6AT1G09200SH3.1 Histone superfamily protein;(source:Araport11)NA
In [48]:
msc <- c()
suppressWarnings(
for (i in as.character(unique(cc.genes.zhang$Phase))){
    if (length(cc.genes.zhang[which(cc.genes.zhang$Phase== i),]$Gene)>1){
    msc <- cbind(msc, as.numeric(apply(apply(as.matrix(rc.integrated@assays$integrated@data)[cc.genes.zhang[which(cc.genes.zhang$Phase== i),]$Gene,], 1, zscore), 1, mean)))       
    } else {
    msc <- cbind(msc, as.numeric(zscore(as.matrix(rc.integrated@assays$integrated@data)[cc.genes.zhang[which(cc.genes.zhang$Phase== i),]$Gene,])))      
    }

}
)
colnames(msc) <- as.character(unique(cc.genes.zhang$Phase))
rownames(msc) <- colnames(rc.integrated)
In [49]:
head(msc)
A matrix: 6 x 3 of type dbl
G2MS
AAACCTGAGCGCCTCA_1-0.028233902 0.02767077-0.2789407
AAACCTGCACTTAACG_1-0.085037611-0.09186232-0.3035347
AAACCTGCAGGACGTA_1-0.128067505-0.03061144-0.2891731
AAACCTGGTATCACCA_1-0.005055825-0.07676504-0.2121654
AAACCTGTCAGGTAAA_1-0.082317857-0.09259665-0.1922454
AAACGGGCAATAGAGT_1-0.154851348-0.11182951-0.3160343
In [51]:
#rc.integrated$G1_sig <- msc[,1]
rc.integrated$G2_sig <- msc[,1]
rc.integrated$M_sig <- msc[,2]
rc.integrated$S_sig <- msc[,3]
In [52]:
#FeaturePlot(rc.integrated, reduction="umap", features = "G1_sig", pt.size=1,order=TRUE)
In [53]:
FeaturePlot(rc.integrated, reduction="umap", features = "S_sig", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [54]:
FeaturePlot(rc.integrated, reduction="umap", features = "G2_sig", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [55]:
FeaturePlot(rc.integrated, reduction="umap", features = "M_sig", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [56]:
#FeaturePlot(rc.integrated, reduction="pca", features = "G1_sig", pt.size=1,order=TRUE)
In [57]:
FeaturePlot(rc.integrated, reduction="pca", features = "S_sig", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [58]:
FeaturePlot(rc.integrated, reduction="pca", features = "G2_sig", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [59]:
FeaturePlot(rc.integrated, reduction="pca", features = "M_sig", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Goldy's list¶

In [60]:
head(cc.genes.goldy)
A data.frame: 6 x 2
genephase
<chr><chr>
1AT2G23320Enriched_G2M
2AT4G17340Enriched_G2M
3AT5G27420Enriched_G2M
4AT5G54290Enriched_G2M
5AT5G08240Enriched_G2M
6AT5G06980Enriched_G2M
In [61]:
cc.genes.goldy <- cc.genes.goldy[match(intersect(rownames(rc.integrated),cc.genes.goldy$gene),cc.genes.goldy$gene),] %>% arrange(phase)
In [62]:
msc <- c()
suppressWarnings(
for (i in as.character(unique(cc.genes.goldy$phase))){
    if (length(cc.genes.goldy[which(cc.genes.goldy$phase== i),]$gene)>1){
    msc <- cbind(msc, as.numeric(apply(apply(as.matrix(rc.integrated@assays$integrated@data)[cc.genes.goldy[which(cc.genes.goldy$phase== i),]$gene,], 1, zscore), 1, mean)))       
    } else {
    msc <- cbind(msc, as.numeric(zscore(as.matrix(rc.integrated@assays$integrated@data)[cc.genes.goldy[which(cc.genes.goldy$phase== i),]$gene,])))      
    }

}
)
colnames(msc) <- as.character(unique(cc.genes.goldy$phase))
rownames(msc) <- colnames(rc.integrated)
In [63]:
head(msc)
A matrix: 6 x 2 of type dbl
Depleted_G2MEnriched_G2M
AAACCTGAGCGCCTCA_1-0.055274931 0.05566898
AAACCTGCACTTAACG_1 0.002564458-0.10651000
AAACCTGCAGGACGTA_1-0.072215184-0.10086421
AAACCTGGTATCACCA_1-0.161443430 0.01979224
AAACCTGTCAGGTAAA_1 0.028097562-0.05966497
AAACGGGCAATAGAGT_1 0.092315593-0.06260356
In [64]:
rc.integrated$Depleted_G2M <- msc[,1]
rc.integrated$Enriched_G2M <- msc[,2]
In [65]:
FeaturePlot(rc.integrated, reduction="umap", features = "Depleted_G2M", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [66]:
FeaturePlot(rc.integrated, reduction="umap", features = "Enriched_G2M", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [67]:
FeaturePlot(rc.integrated, reduction="pca", features = "Depleted_G2M", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [68]:
FeaturePlot(rc.integrated, reduction="pca", features = "Enriched_G2M", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Check Seurat Cluster 0 and 4¶

In [76]:
rc.integrated$G2M_clust <- rep("others", ncol(rc.integrated))
rc.integrated$G2M_clust[which(rc.integrated$seurat_clusters==4)] = "G2M"
In [77]:
DimPlot(rc.integrated, reduction="umap", group = "G2M_clust", pt.size=1,order=FALSE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [81]:
DimPlot(rc.integrated, reduction="pca", group = "G2M_clust", pt.size=1,order=FALSE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [78]:
rc.integrated$S_clust <- rep("others", ncol(rc.integrated))
rc.integrated$S_clust[which(rc.integrated$seurat_clusters==0)] = "S"
In [79]:
DimPlot(rc.integrated, reduction="umap", group = "S_clust", pt.size=1,order=FALSE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [82]:
DimPlot(rc.integrated, reduction="pca", group = "S_clust", pt.size=1,order=FALSE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Convert to single cell experiment¶

In [10]:
new.ref <- readRDS("./tradeseq/Tricycle_Reference_Arabidopsis_Root.rds")
In [11]:
sce <- as.SingleCellExperiment(rc.integrated)
In [12]:
sce$geno <- rep("WT",length(sce$orig.ident))
In [13]:
## Tricycle estimation
sce <- project_cycle_space(sce,ref.m = new.ref)
sce <- estimate_cycle_position(sce)

sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
sce$tricyclePosition <- sce$tricyclePosition*pi

## Plot
options(repr.plot.width=9, repr.plot.height=4)
print(plot_ccposition_den(sce$tricyclePosition,
                sce$geno, 'genotype',
                bw = 10, fig.title = "Kernel density of \u03b8", palette.v=colorRampPalette(brewer.pal(11, "Spectral"))(18)) +
  theme_bw(base_size = 14)+ geom_vline(xintercept = 1*pi, linetype="dotted", 
                color = "red", size=1.5)+ geom_vline(xintercept = 1.55*pi, linetype="dotted", 
                color = "red", size=1.5))
The number of projection genes found in the new data is 245.

Warning message:
“Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.”
In [14]:
grDevices::cairo_pdf("./pdfs/Density_Plot_CC.pdf",width=9, height=4)
print(plot_ccposition_den(sce$tricyclePosition,
                sce$geno, 'genotype',
                bw = 10, fig.title = "Kernel density of \u03b8", palette.v=colorRampPalette(brewer.pal(11, "Spectral"))(18)) +
  theme_bw(base_size = 14)+ geom_vline(xintercept = 1*pi, linetype="dotted", 
                color = "red", size=1.5)+ geom_vline(xintercept = 1.55*pi, linetype="dotted", 
                color = "red", size=1.5))
dev.off()
png: 2
In [16]:
p <- plot_emb_circle_scale(sce, dimred = 1,
                           point.size = 5, point.alpha = 0.9) +
  theme_bw(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
options(repr.plot.width=8, repr.plot.height=4)
plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.6))
In [17]:
grDevices::cairo_pdf("./pdfs/PCA_CC.pdf",width=8, height=4)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.6)))
dev.off()
png: 2
In [16]:
## Tricycle estimation
sce <- project_cycle_space(sce,ref.m = new.ref)
sce <- estimate_cycle_position(sce)

sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
sce$tricyclePosition <- sce$tricyclePosition*pi

## Plot
options(repr.plot.width=9, repr.plot.height=9)
print(plot_ccposition_den(sce$tricyclePosition,
                sce$orig.ident, 'sample',
                bw = 10, fig.title = "Kernel density of \u03b8", palette.v=colorRampPalette(brewer.pal(11, "Spectral"))(18)) +
  theme_bw(base_size = 14)+ geom_vline(xintercept = 1*pi, linetype="dotted", 
                color = "red", size=1.5)+ geom_vline(xintercept = 1.55*pi, linetype="dotted", 
                color = "red", size=1.5))
The number of projection genes found in the new data is 245.

In [85]:
## Assign position to stage
    tricycleCCStage <- sce$tricyclePosition/pi

    G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
    S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
    G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
                    
    tricycleCCStage[G1_idx]="G1/G0"
    tricycleCCStage[S_idx]="S"
    tricycleCCStage[G2_idx]="G2M"
    
    sce$tricycleCCStage <- tricycleCCStage
    rc.integrated$tricyclePosition <- sce$tricyclePosition
    rc.integrated$tricycleCCStage <- tricycleCCStage
In [92]:
saveRDS(rc.integrated,"./scRNA-seq/Integrated_Objects/rc.integrated_18S_WT_cell_cycle_seu4_245genes_Goldy_Zhang_AMIGO_20230825.rds")
In [86]:
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                       point.size = 1, point.alpha = 0.9) +
  theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
In [87]:
DimPlot(rc.integrated, reduction="umap", group = "tricycleCCStage", pt.size=1,order=FALSE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [88]:
DimPlot(rc.integrated, reduction="pca", group = "tricycleCCStage", pt.size=1,order=FALSE)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

In [89]:
rc.integrated$CCclust <- rep("others", ncol(rc.integrated))
rc.integrated$CCclust[which(rc.integrated$seurat_clusters==0)] = "S"
rc.integrated$CCclust[which(rc.integrated$seurat_clusters==4)] = "G2M"
table(rc.integrated$tricycleCCStage, rc.integrated$CCclust)
       
           G2M others      S
  G1/G0    920 154117   1510
  G2M     6101   4126    117
  S        143  20007  12769
In [90]:
scater::plotReducedDim(sce, dimred = "tricycleEmbedding") +
    labs(x = "Projected PC1", y = "Projected PC2") +
    ggtitle(sprintf("Projected cell cycle space (n=%d)",
                    ncol(sce))) +
    theme_bw(base_size = 14)
In [91]:
p <- plot_emb_circle_scale(sce, dimred = "tricycleEmbedding",
                           point.size = 3.5, point.alpha = 0.9) +
  theme_bw(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.4))
In [ ]:
## Total RNA content
#fit.l <- fit_periodic_loess(sce$tricyclePosition,
#                            log2(sce$nCount_RNA),
#                            plot = TRUE, span=0.3,
#                       x_lab = "Cell cycle position \u03b8", y_lab = "log2(Total UMIs)",
#                       fig.title = paste0("Expression of key gene along \u03b8 (n=",
#                                          ncol(sce), ")"))

#fit <- fit.l$pred.df
#options(repr.plot.width=8, repr.plot.height=6)
#ggplot(fit, aes(x=x,y=y)) + geom_path(size = 1, alpha = 1)+ labs(x = "θ", y = "log2(total UMIs)", title = "") + 
#            scale_x_continuous(limits = c(0, 2 * pi), breaks = c(0, 
#                pi/2, pi, 3 * pi/2, 2 * pi), labels = paste0(c(0, 
#                0.5, 1, 1.5, 2), "π"), name = "θ")+
#  theme_bw(base_size = 14)+
#  theme(legend.title=element_blank())
In [11]:
table(rc.integrated$time.anno.Li)
         Distal Columella   Distal Lateral Root Cap                Elongation 
                     6526                      4904                     83690 
               Maturation      Proliferation Domain        Proximal Columella 
                    21054                     21913                      3063 
Proximal Lateral Root Cap         Transition Domain 
                    27495                     31165 
In [12]:
sub <- subset(rc.integrated, cells=colnames(rc.integrated)[which(rc.integrated$time.anno.Li=="Proliferation Domain"|
                                                                rc.integrated$time.anno.Li=="Proximal Columella"|
                                                                rc.integrated$time.anno.Li=="Proximal Lateral Root Cap")])
In [13]:
sub
An object of class Seurat 
57118 features across 52471 samples within 3 assays 
Active assay: integrated (245 features, 245 variable features)
 2 layers present: data, scale.data
 2 other assays present: RNA, SCT
 3 dimensional reductions calculated: pca, umap, umap_2D
In [14]:
table(sub$time.anno.Li)
     Proliferation Domain        Proximal Columella Proximal Lateral Root Cap 
                    21913                      3063                     27495 
In [15]:
# Run PCA
sub <- RunPCA(sub, npcs = 10, verbose = FALSE, approx = TRUE)

# Find nearest neighbors
#rc.integrated <- FindNeighbors(rc.integrated, reduction = "pca",dims = 1:50)
Warning message:
“Number of dimensions changing from 50 to 10”
Warning message:
“The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session”
20:27:45 UMAP embedding parameters a = 0.9922 b = 1.112

Found more than one class "dist" in cache; using the first, from namespace 'spam'

Also defined by ‘BiocGenerics’

20:27:45 Read 52471 rows and found 5 numeric columns

20:27:45 Using Annoy for neighbor search, n_neighbors = 30

Found more than one class "dist" in cache; using the first, from namespace 'spam'

Also defined by ‘BiocGenerics’

20:27:45 Building Annoy index with metric = cosine, n_trees = 50

0%   10   20   30   40   50   60   70   80   90   100%

[----|----|----|----|----|----|----|----|----|----|

*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
|

20:27:57 Writing NN index file to temp file /tmp/Rtmpu0VjyG/file317db91a27d9c1

20:27:57 Searching Annoy index using 1 thread, search_k = 3000

20:28:40 Annoy recall = 100%

20:28:42 Commencing smooth kNN distance calibration using 1 thread
 with target n_neighbors = 30

20:28:44 5 smooth knn distance failures

20:28:51 Initializing from normalized Laplacian + noise (using RSpectra)

20:28:59 Commencing optimization for 200 epochs, with 1889004 positive edges

20:29:42 Optimization finished

Warning message:
“Key ‘UMAP_’ taken, using ‘umap_’ instead”
In [23]:
# Run 2D UMAP
sub <- RunUMAP(sub, reduction = "pca", dims = 1:3, n.neighbors = 30, min.dist=0.3)
20:40:15 UMAP embedding parameters a = 0.9922 b = 1.112

Found more than one class "dist" in cache; using the first, from namespace 'spam'

Also defined by ‘BiocGenerics’

20:40:15 Read 52471 rows and found 3 numeric columns

20:40:15 Using Annoy for neighbor search, n_neighbors = 30

Found more than one class "dist" in cache; using the first, from namespace 'spam'

Also defined by ‘BiocGenerics’

20:40:15 Building Annoy index with metric = cosine, n_trees = 50

0%   10   20   30   40   50   60   70   80   90   100%

[----|----|----|----|----|----|----|----|----|----|

*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
|

20:40:28 Writing NN index file to temp file /tmp/Rtmpu0VjyG/file317db92623cbe

20:40:28 Searching Annoy index using 1 thread, search_k = 3000

20:40:57 Annoy recall = 100%

20:40:59 Commencing smooth kNN distance calibration using 1 thread
 with target n_neighbors = 30

20:41:00 18 smooth knn distance failures

20:41:06 Initializing from normalized Laplacian + noise (using RSpectra)

20:41:40 Commencing optimization for 200 epochs, with 1402752 positive edges

20:42:13 Optimization finished

Warning message:
“Key ‘UMAP_’ taken, using ‘umap_’ instead”
In [49]:
my_cols = brewer.pal(8,"Set2")
In [50]:
options(repr.plot.width=10, repr.plot.height=8)
DimPlot(sub, reduction = "umap", group.by = "time.anno.Li",cols=alpha(my_cols,0.8), label = FALSE, pt.size=1)
In [24]:
sce <- as.SingleCellExperiment(sub)
In [30]:
# dims = 3
options(repr.plot.width=8, repr.plot.height=8)
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                       point.size = 2, point.alpha = 0.8) +
  theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 1)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
In [22]:
# dims = 4
options(repr.plot.width=8, repr.plot.height=8)
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                       point.size = 2, point.alpha = 1) +
  theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 1)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
In [19]:
# dims = 5
options(repr.plot.width=8, repr.plot.height=8)
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                       point.size = 2, point.alpha = 1) +
  theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 1)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
In [ ]:
# dims = 5
options(repr.plot.width=8, repr.plot.height=8)
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                       point.size = 2, point.alpha = 1) +
  theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 1)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))

Tricycle inference¶

In [91]:
run_tricycle <- function(sample){
    seu <- readRDS(paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT.rds"))
    sce <- as.SingleCellExperiment(seu)
    
    ## Tricycle estimation
    sce <- project_cycle_space(sce,ref.m = new.ref)
    sce <- estimate_cycle_position(sce)
    
    sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
    sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
    sce$tricyclePosition <- sce$tricyclePosition*pi
    
    ## Plot
    print(plot_ccposition_den(sce$tricyclePosition,
                    sce$orig.ident, 'sample',
                    bw = 10, fig.title = "Kernel density of \u03b8") +
      theme_bw(base_size = 14))
    
    ## Assign position to stage
    tricycleCCStage <- sce$tricyclePosition/pi

    G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
    S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
    G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
                    
    tricycleCCStage[G1_idx]="G1/G0"
    tricycleCCStage[S_idx]="S"
    tricycleCCStage[G2_idx]="G2M"
    
    sce$tricycleCCStage <- tricycleCCStage
    
    ## Schwabe Stage
    #sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
    
    ## Table
    #seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
    #print(table(seu$time.anno.Li, sce$tricycleCCStage))
    
    ## Prepare seu
    #seu <- prep(seu)
    #plot_anno(seu)
    
    ## Sce to seu
    seu$tricyclePosition <- sce$tricyclePosition/pi
    seu$tricycleCCStage <- sce$tricycleCCStage
    #seu$SchwabeCCStage <- sce$CCStage
    
    ## Tricycle position on umap
    p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                           point.size = 1, point.alpha = 0.9) +
      theme_void(base_size = 14)
    legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
    print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
    
    ## Write annotations as csv
    write.csv(seu@meta.data[,c("tricyclePosition", "tricycleCCStage")], paste0("./tricycle/",sample,".csv"), quote=FALSE)
    #saveRDS(seu,paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT_seu4.rds"))
}
In [92]:
sample.list <- c("sc_130","sc_131","sc_132","sc_133","sc_134","sc_135","sc_136","sc_137")
In [93]:
for (i in 1:length(sample.list)){
 run_tricycle(sample.list[i])
}
The number of projection genes found in the new data is 230.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 237.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
In [94]:
sample.list <- c("sc_5")
In [95]:
for (i in 1:length(sample.list)){
 run_tricycle(sample.list[i])
}
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
In [96]:
sample.list <- c("sc_1","sc_2","sc_111","sc_112","sc_113","sc_114")
In [97]:
for (i in 1:length(sample.list)){
 run_tricycle(sample.list[i])
}
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 243.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
In [98]:
sample.list <- c("sc_151","sc_152","sc_153","sc_154")
In [99]:
for (i in 1:length(sample.list)){
 run_tricycle(sample.list[i])
}
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
In [100]:
sample.list <- c("pp1","col0","tnw2","sc_1","sc_9_at","sc_30", "sc_31", "sc_37", "sc_40", "sc_51")
In [101]:
for (i in 1:length(sample.list)){
 run_tricycle(sample.list[i])
}
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 243.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 244.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 241.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
In [ ]:
run_tricycle <- function(sample){
    seu <- readRDS(paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT.rds"))
    sce <- as.SingleCellExperiment(seu)
    
    ## Tricycle estimation
    sce <- project_cycle_space(sce,ref.m = new.ref)
    sce <- estimate_cycle_position(sce)
    
    sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.2
    sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
    sce$tricyclePosition <- (sce$tricyclePosition-2)*-1
    sce$tricyclePosition <- sce$tricyclePosition*pi
    
    ## Plot
    print(plot_ccposition_den(sce$tricyclePosition,
                    sce$orig.ident, 'sample',
                    bw = 10, fig.title = "Kernel density of \u03b8") +
      theme_bw(base_size = 14))
    
    ## Assign position to stage
    tricycleCCStage <- sce$tricyclePosition/pi

    G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
    S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
    G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
                    
    tricycleCCStage[G1_idx]="G1/G0"
    tricycleCCStage[S_idx]="S"
    tricycleCCStage[G2_idx]="G2M"
    
    sce$tricycleCCStage <- tricycleCCStage
    
    ## Schwabe Stage
    #sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
    
    ## Table
    #seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
    #print(table(seu$time.anno.Li, sce$tricycleCCStage))
    
    ## Prepare seu
    #seu <- prep(seu)
    #plot_anno(seu)
    
    ## Sce to seu
    seu$tricyclePosition <- sce$tricyclePosition/pi
    seu$tricycleCCStage <- sce$tricycleCCStage
    #seu$SchwabeCCStage <- sce$CCStage
    
    ## Tricycle position on umap
    p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                           point.size = 1, point.alpha = 0.9) +
      theme_void(base_size = 14)
    legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
    print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
    
    ## Write annotations as csv
    write.csv(seu@meta.data[,c("tricyclePosition", "tricycleCCStage")], paste0("./tricycle/",sample,".csv"), quote=FALSE)
    #saveRDS(seu,paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT_seu4.rds"))
}
In [94]:
sample.list <- c("dc1","dc2","tnw1","briT","briTR","sc_10_at","sc_11","sc_12","sc_20","sc_43","sc_44","sc_45","sc_46","sc_47","sc_48","sc_49","sc_50","sc_71",
                 "sc_122","sc_123","sc_124","sc_125","sc_126","sc_127","sc_128","sc_129","sc_130","sc_131","sc_132","sc_133","sc_134","sc_135","sc_136","sc_137",
                 "sc_155","sc_156","sc_157","sc_158","sc_159","sc_160","sc_161","sc_162","sc_163","sc_164","sc_165","sc_166","sc_167","sc_168","sc_169",
                 "sc_170","sc_171","sc_172","sc_173","sc_175","sc_176","sc_177","sc_178","sc_179","sc_180","sc_186","sc_187","sc_189","sc_190","sc_191",
                 "sc_215","sc_216","sc_217","sc_218","sc_220","sc_221","sc_222","sc_223","sc_224","sc_225","sc_226","sc_227","sc_228","sc_229","sc_235")
In [95]:
for (i in 1:length(sample.list)){
 run_tricycle(sample.list[i])
}
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 240.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 228.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 244.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

Permanent recipe for tricycle inference:¶

In [21]:
seu <- readRDS(paste0("./scRNA-seq/Seurat_Objects/","dc1","_COPILOT.rds"))
In [22]:
sce <- as.SingleCellExperiment(seu)
    
    ## Tricycle estimation
    sce <- project_cycle_space(sce,ref.m = new.ref)
    sce <- estimate_cycle_position(sce)
The number of projection genes found in the new data is 245.

In [23]:
sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
    sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
    sce$tricyclePosition <- sce$tricyclePosition*pi
    
    ## Plot
    print(plot_ccposition_den(sce$tricyclePosition,
                    sce$orig.ident, 'sample',
                    bw = 10, fig.title = "Kernel density of \u03b8") +
      theme_bw(base_size = 14))
    
    ## Assign position to stage
    tricycleCCStage <- sce$tricyclePosition/pi

    G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
    S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
    G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
                    
    tricycleCCStage[G1_idx]="G1/G0"
    tricycleCCStage[S_idx]="S"
    tricycleCCStage[G2_idx]="G2M"
    
    sce$tricycleCCStage <- tricycleCCStage
    
    ## Schwabe Stage
    #sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
    
    ## Table
    #seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
    #print(table(seu$time.anno.Li, sce$tricycleCCStage))
    
    ## Prepare seu
    #seu <- prep(seu)
    #plot_anno(seu)
    
    ## Sce to seu
    seu$tricyclePosition <- sce$tricyclePosition/pi
    seu$tricycleCCStage <- sce$tricycleCCStage
    #seu$SchwabeCCStage <- sce$CCStage
    
    ## Tricycle position on umap
    p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                           point.size = 1, point.alpha = 0.9) +
      theme_void(base_size = 14)
    legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
    print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
In [24]:
run_tricycle <- function(sample){
    seu <- readRDS(paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT.rds"))
    sce <- as.SingleCellExperiment(seu)
    
    ## Tricycle estimation
    sce <- project_cycle_space(sce,ref.m = new.ref)
    sce <- estimate_cycle_position(sce)
    
    sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
    sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
    sce$tricyclePosition <- sce$tricyclePosition*pi
    
    ## Plot
    print(plot_ccposition_den(sce$tricyclePosition,
                    sce$orig.ident, 'sample',
                    bw = 10, fig.title = "Kernel density of \u03b8") +
      theme_bw(base_size = 14))
    
    ## Assign position to stage
    tricycleCCStage <- sce$tricyclePosition/pi

    G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
    S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
    G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
                    
    tricycleCCStage[G1_idx]="G1/G0"
    tricycleCCStage[S_idx]="S"
    tricycleCCStage[G2_idx]="G2M"
    
    sce$tricycleCCStage <- tricycleCCStage
    
    ## Schwabe Stage
    #sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
    
    ## Table
    #seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
    #print(table(seu$time.anno.Li, sce$tricycleCCStage))
    
    ## Prepare seu
    #seu <- prep(seu)
    #plot_anno(seu)
    
    ## Sce to seu
    seu$tricyclePosition <- sce$tricyclePosition/pi
    seu$tricycleCCStage <- sce$tricycleCCStage
    #seu$SchwabeCCStage <- sce$CCStage
    
    ## Tricycle position on umap
    p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                           point.size = 1, point.alpha = 0.9) +
      theme_void(base_size = 14)
    legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
    print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
    
    ## Write annotations as csv
    write.csv(seu@meta.data[,c("tricyclePosition", "tricycleCCStage")], paste0("./tricycle/",sample,".csv"), quote=FALSE)
    #saveRDS(seu,paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT_seu4.rds"))
}

Protophloem SMART-seq data¶

In [25]:
mtx <- read.table('./scRNA-seq/Published_Data/GSE140778_20190115_data.txt')
In [26]:
meta <- read.csv('./scRNA-seq/Published_Data/GSE140778_20190115_meta_data.csv')
In [27]:
sum(!is.na(meta$Cluster_PSE))
758
In [28]:
head(mtx)
A data.frame: 6 x 1375
SLX_12100_i701_i502SLX_12100_i701_i503SLX_12100_i701_i504SLX_12100_i701_i505SLX_12100_i701_i506SLX_12100_i701_i507SLX_12100_i701_i508SLX_12100_i701_i517SLX_12100_i702_i502SLX_12100_i702_i503...SLX_17313_i728_i521SLX_17313_i728_i522SLX_17313_i729_i513SLX_17313_i729_i515SLX_17313_i729_i516SLX_17313_i729_i517SLX_17313_i729_i518SLX_17313_i729_i520SLX_17313_i729_i521SLX_17313_i729_i522
<dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>...<dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
AT1G01010138.832 0.0000 1.37359 0.665616 23.122100 0.000000.00000153.1940...45.11770 0.000000 13.754358.2903 0.000.000000 0.00000
AT1G01020277.664146.593033.42390 0.000000 0.000000 39.605400.00000105.3210... 8.01155 0.000000 0.0000 0.0000 0.000.000000 0.00000
AT1G01030 0.000 0.0000 0.00000 7.321770 0.000000 0.000000.00000 0.0000... 0.00000 0.000000 0.0000 0.0000 0.000.000000 0.00000
AT1G01040 0.000 0.0000 0.00000 0.332834 0.379051 22.996700.00000 4.7873... 0.00000 2.438140 0.000026.1001 0.000.000000 0.00000
AT1G01046 0.000 0.0000 0.00000 0.000000 0.000000 0.000000.00000 0.0000... 0.00000 0.000000 0.0000 0.0000 0.000.000000 0.00000
AT1G01050138.832146.593072.80000105.167000137.974000136.702001.94657 4.7873... 6.32491104.809000103.157029.5801329.676.51429010.54810
In [29]:
rownames(meta) <- meta$ID
In [30]:
seu <- CreateSeuratObject(mtx, meta.data=meta)
In [31]:
seu
An object of class Seurat 
32824 features across 1375 samples within 1 assay 
Active assay: RNA (32824 features, 0 variable features)
In [32]:
head(seu@meta.data)
A data.frame: 6 x 16
orig.identnCount_RNAnFeature_RNAIDsamplereads_totalreads_unmappedreads_uniquereads_multipercent_uniquepercent_multipercent_mappedlib_size_avrCluster_all_cellsCluster_PSEPseudotime_PSE
<fct><dbl><int><chr><chr><int><int><int><int><dbl><dbl><dbl><int><int><int><dbl>
SLX_12100_i701_i502SLX734281.9 2860SLX_12100_i701_i502P1 847809 803895 38147 5767 4.499480 0.6802240 5.179704678NANA NA
SLX_12100_i701_i503SLX913422.6 3518SLX_12100_i701_i503P114316591399319 27698 4642 1.934679 0.3242392 2.258918678NANA NA
SLX_12100_i701_i504SLX827476.9 2904SLX_12100_i701_i504P1 313281 300049 9440 3792 3.013269 1.2104149 4.223684678NANA NA
SLX_12100_i701_i505SLX600775.510725SLX_12100_i701_i505P13328241 9923321515355 82055445.53020624.654284470.18449167810 41.3333992
SLX_12100_i701_i506SLX922948.511363SLX_12100_i701_i506P1705164319010223201903194871845.40648227.634949873.04143267814 40.8876946
SLX_12100_i701_i507SLX894200.010634SLX_12100_i701_i507P1690631018814992873579215123241.60802231.148790072.756812678 7NA NA
In [33]:
seu <- NormalizeData(seu)
seu <- FindVariableFeatures(seu, selection.method = "vst", nfeatures = nrow(seu))
all.genes <- rownames(seu)
seu <- ScaleData(seu, features = all.genes)
seu <- RunPCA(seu,features = VariableFeatures(object = seu))
Centering and scaling data matrix

PC_ 1 
Positive:  AT5G03240, AT3G53990, AT4G05320, AT1G70670, AT4G05050, AT1G75280, AT4G39090, AT3G22400, AT4G21980, AT1G29520 
	   AT3G03200, AT1G77700, AT5G48060, AT1G54400, AT3G06420, AT3G58720, AT4G33490, AT2G46630, AT1G29380, AT1G08160 
	   AT1G79480, AT5G07220, AT5G49320, AT5G18970, AT5G03510, AT2G41200, AT2G42820, AT1G76970, AT1G63310, AT1G61760 
Negative:  AT3G04400, AT3G24830, AT3G25520, AT5G20290, AT5G02960, AT1G08360, AT5G16130, AT1G18540, AT3G02080, AT3G49910 
	   AT4G27090, AT2G34480, AT5G15200, AT2G47610, AT5G10360, AT4G00100, AT4G34670, AT5G23740, AT2G39460, AT3G04920 
	   AT3G53020, AT1G69620, AT2G18020, AT2G27530, AT1G67430, AT3G14600, AT3G11940, AT3G09200, AT2G19730, AT4G18100 
PC_ 2 
Positive:  AT5G42050, AT1G20450, AT1G19180, AT1G20440, AT2G47730, AT1G78380, AT4G30960, AT3G10985, AT5G10830, AT5G11670 
	   AT5G45110, AT1G78660, AT5G63790, AT1G04250, AT2G29440, AT5G08790, AT4G16380, AT4G34150, AT5G64260, AT5G19110 
	   AT2G24850, AT4G35090, AT4G19200, AT3G52470, AT5G06320, AT4G21990, AT1G55920, AT1G78080, AT1G76680, AT4G36040 
Negative:  AT5G62960, AT1G80940, AT2G38550, AT1G05610, AT2G22730, AT1G31880, AT5G03530, AT1G55760, AT2G18380, AT1G62045 
	   AT1G79430, AT5G40580, AT1G63310, AT1G68000, AT1G61760, AT3G49360, AT4G23900, AT1G69970, AT2G20515, AT2G36400 
	   AT4G17510, AT5G26940, AT1G14900, AT5G04060, AT1G06490, AT2G32280, AT5G17260, AT5G48060, AT4G24250, AT1G79480 
PC_ 3 
Positive:  AT1G03220, AT1G48920, AT1G66580, AT5G52470, AT3G15000, AT5G20160, AT3G10690, AT1G74560, AT4G12600, AT1G56110 
	   AT5G27120, AT4G26110, AT2G38810, AT5G25370, AT4G25630, AT1G68230, AT3G05060, AT3G07770, AT3G44750, AT5G22650 
	   AT3G45970, AT1G52930, AT3G23990, AT3G58350, AT3G56070, AT2G19540, AT3G57150, AT4G15910, AT3G03950, AT1G26680 
Negative:  AT1G78040, AT5G60920, AT5G59880, AT2G45470, AT3G12610, AT5G39320, AT5G05170, AT5G49720, AT1G53500, AT1G75680 
	   AT1G08200, AT4G12730, AT5G12250, AT2G20750, AT1G26850, AT1G70490, AT1G19835, AT4G31590, AT3G12710, AT1G10630 
	   AT1G76160, AT4G34490, AT1G45688, AT4G32410, AT1G76670, AT5G59290, AT5G64740, AT5G47770, AT5G25100, AT3G07010 
PC_ 4 
Positive:  AT1G10760, AT4G01610, AT1G54030, AT1G32860, AT2G34020, AT1G18210, AT2G29980, AT3G55760, AT5G48300, AT3G12490 
	   AT4G09020, AT4G31290, AT1G21340, AT3G08690, AT1G56680, AT2G15490, AT2G32430, AT2G04025, AT4G00490, AT5G08590 
	   AT2G28900, AT3G46440, AT5G18470, AT2G41800, AT4G34970, AT1G72490, AT1G13620, AT4G23060, AT3G18230, AT1G19850 
Negative:  AT4G14130, AT1G62480, AT4G21960, AT3G53420, AT2G02130, AT5G59090, AT1G52760, AT3G21770, AT2G43910, AT3G04730 
	   AT4G29100, AT3G19450, AT3G23820, AT3G51950, AT5G56540, AT3G14310, AT4G11210, AT3G56240, AT2G30520, AT3G16330 
	   AT5G24800, AT4G19410, AT2G31083, AT5G51780, AT2G44790, AT1G10682, AT3G16340, AT3G11930, AT5G25890, AT4G36710 
PC_ 5 
Positive:  AT4G15960, AT1G70670, AT3G02470, AT4G18710, AT2G32150, AT2G39420, AT2G36650, AT3G14780, AT4G37730, AT5G16650 
	   AT3G47160, AT4G35450, AT4G35790, AT2G38800, AT5G13740, AT3G11280, AT4G39090, AT2G19110, AT5G64572, AT1G75380 
	   AT5G49320, AT1G80040, AT1G44170, AT5G18980, AT3G60680, AT3G63000, AT3G14350, AT4G21110, AT2G39080, AT2G21270 
Negative:  ATMG01390, ATMG00020, AT5G18610, AT2G07698, AT1G77690, AT5G23860, ATMG01190, AT1G04680, AT5G26270, AT3G02230 
	   AT2G04780, ATMG01360, AT1G64640, AT3G54810, AT1G69690, AT2G14890, ATCG00500, AT2G20750, ATMG00660, AT2G47650 
	   ATCG01180, ATCG00470, AT3G48185, AT2G21160, AT3G58120, ATCG00950, AT5G15650, ATMG00090, AT2G34300, AT1G21750 

In [34]:
seu <- FindNeighbors(seu, dims = 1:10)
seu <- FindClusters(seu, resolution = 0.5)
Computing nearest neighbor graph

Computing SNN

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 1375
Number of edges: 33789

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9051
Number of communities: 12
Elapsed time: 0 seconds
In [35]:
seu <- RunUMAP(seu, dims = 1:10)
Warning message:
“The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session”
10:49:41 UMAP embedding parameters a = 0.9922 b = 1.112

10:49:41 Read 1375 rows and found 10 numeric columns

10:49:41 Using Annoy for neighbor search, n_neighbors = 30

10:49:41 Building Annoy index with metric = cosine, n_trees = 50

0%   10   20   30   40   50   60   70   80   90   100%

[----|----|----|----|----|----|----|----|----|----|

*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
|

10:49:41 Writing NN index file to temp file /tmp/RtmpSzDebc/file16ca8b2394d636

10:49:41 Searching Annoy index using 1 thread, search_k = 3000

10:49:42 Annoy recall = 100%

10:49:43 Commencing smooth kNN distance calibration using 1 thread
 with target n_neighbors = 30

10:49:46 Initializing from normalized Laplacian + noise (using RSpectra)

10:49:46 Commencing optimization for 500 epochs, with 49356 positive edges

10:49:49 Optimization finished

In [36]:
# note that you can set `label = TRUE` or use the LabelClusters function to help label
# individual clusters
DimPlot(seu, reduction = "umap")
In [37]:
head(seu@meta.data)
A data.frame: 6 x 18
orig.identnCount_RNAnFeature_RNAIDsamplereads_totalreads_unmappedreads_uniquereads_multipercent_uniquepercent_multipercent_mappedlib_size_avrCluster_all_cellsCluster_PSEPseudotime_PSERNA_snn_res.0.5seurat_clusters
<fct><dbl><int><chr><chr><int><int><int><int><dbl><dbl><dbl><int><int><int><dbl><fct><fct>
SLX_12100_i701_i502SLX734281.9 2860SLX_12100_i701_i502P1 847809 803895 38147 5767 4.499480 0.6802240 5.179704678NANA NA1010
SLX_12100_i701_i503SLX913422.6 3518SLX_12100_i701_i503P114316591399319 27698 4642 1.934679 0.3242392 2.258918678NANA NA4 4
SLX_12100_i701_i504SLX827476.9 2904SLX_12100_i701_i504P1 313281 300049 9440 3792 3.013269 1.2104149 4.223684678NANA NA0 0
SLX_12100_i701_i505SLX600775.510725SLX_12100_i701_i505P13328241 9923321515355 82055445.53020624.654284470.18449167810 41.33339920 0
SLX_12100_i701_i506SLX922948.511363SLX_12100_i701_i506P1705164319010223201903194871845.40648227.634949873.04143267814 40.88769468 8
SLX_12100_i701_i507SLX894200.010634SLX_12100_i701_i507P1690631018814992873579215123241.60802231.148790072.756812678 7NA NA7 7
In [38]:
seu
An object of class Seurat 
32824 features across 1375 samples within 1 assay 
Active assay: RNA (32824 features, 27218 variable features)
 2 dimensional reductions calculated: pca, umap
In [39]:
#seu[["umap"]]@cell.embeddings[,1] <- seu[["umap"]]@cell.embeddings[,1]*-1
seu[["umap"]]@cell.embeddings[,2] <- seu[["umap"]]@cell.embeddings[,2]*-1
#u2 <- rc.integrated@reductions$umap@cell.embeddings[,1]
#u1 <- rc.integrated@reductions$umap@cell.embeddings[,2]
#rc.integrated@reductions$umap@cell.embeddings[,1] <- u1
#rc.integrated@reductions$umap@cell.embeddings[,2] <- u2
In [40]:
DimPlot(seu, reduction = "umap", group.by = "Cluster_PSE")
In [41]:
FeaturePlot(seu, reduction = "umap", features = "Pseudotime_PSE")

tricycle reference¶

In [42]:
sce <- as.SingleCellExperiment(seu)

## Tricycle estimation
sce <- project_cycle_space(sce,ref.m = new.ref)
sce <- estimate_cycle_position(sce)
The number of projection genes found in the new data is 245.

In [43]:
sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
    sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
    sce$tricyclePosition <- sce$tricyclePosition*pi
    
    ## Plot
    print(plot_ccposition_den(sce$tricyclePosition,
                    sce$orig.ident, 'sample',
                    bw = 10, fig.title = "Kernel density of \u03b8") +
      theme_bw(base_size = 14))
    
    ## Assign position to stage
    tricycleCCStage <- sce$tricyclePosition/pi

    G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
    S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
    G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
                    
    tricycleCCStage[G1_idx]="G1/G0"
    tricycleCCStage[S_idx]="S"
    tricycleCCStage[G2_idx]="G2M"
    
    sce$tricycleCCStage <- tricycleCCStage
    
    ## Schwabe Stage
    #sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
    
    ## Table
    #seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
    #print(table(seu$time.anno.Li, sce$tricycleCCStage))
    
    ## Prepare seu
    #seu <- prep(seu)
    #plot_anno(seu)
    
    ## Sce to seu
    seu$tricyclePosition <- sce$tricyclePosition/pi
    seu$tricycleCCStage <- sce$tricycleCCStage
    #seu$SchwabeCCStage <- sce$CCStage
Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
In [44]:
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                       point.size = 4, point.alpha = 1) +
  theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 1)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
In [45]:
DimPlot(seu, reduction = "umap", group.by = "tricycleCCStage")
In [46]:
seu$Pseudotime_PSE
SLX_12100_i701_i502
<NA>
SLX_12100_i701_i503
<NA>
SLX_12100_i701_i504
<NA>
SLX_12100_i701_i505
1.333399152
SLX_12100_i701_i506
0.887694607
SLX_12100_i701_i507
<NA>
SLX_12100_i701_i508
1.162116692
SLX_12100_i701_i517
<NA>
SLX_12100_i702_i502
<NA>
SLX_12100_i702_i503
<NA>
SLX_12100_i702_i504
0.092332697
SLX_12100_i702_i505
<NA>
SLX_12100_i702_i506
0.582517359
SLX_12100_i702_i507
2.337060737
SLX_12100_i702_i508
0.854771632
SLX_12100_i702_i517
<NA>
SLX_12100_i703_i502
<NA>
SLX_12100_i703_i503
<NA>
SLX_12100_i703_i504
<NA>
SLX_12100_i703_i505
53.61314228
SLX_12100_i703_i506
54.52189408
SLX_12100_i703_i507
85.19729857
SLX_12100_i703_i508
1.210343071
SLX_12100_i703_i517
14.52197083
SLX_12100_i704_i502
13.24045026
SLX_12100_i704_i503
1.778844138
SLX_12100_i704_i504
52.82139272
SLX_12100_i704_i505
0.437802713
SLX_12100_i704_i506
<NA>
SLX_12100_i704_i507
<NA>
SLX_12100_i704_i508
9.617037557
SLX_12100_i704_i517
52.8213696
SLX_12100_i705_i502
27.10965105
SLX_12100_i705_i503
49.37018009
SLX_12100_i705_i504
1.153673409
SLX_12100_i705_i505
39.32187208
SLX_12100_i705_i506
<NA>
SLX_12100_i705_i507
40.04760222
SLX_12100_i705_i508
26.91993546
SLX_12100_i705_i517
49.9319981
SLX_12100_i706_i502
52.82134648
SLX_12100_i706_i503
27.64106385
SLX_12100_i706_i504
0.626094315
SLX_12100_i706_i506
<NA>
SLX_12100_i706_i507
53.11332564
SLX_12100_i706_i508
<NA>
SLX_12100_i706_i517
49.68501244
SLX_12100_i707_i502
6.766869249
SLX_12100_i707_i503
<NA>
SLX_12100_i707_i504
<NA>
SLX_12100_i707_i505
<NA>
SLX_12100_i707_i506
<NA>
SLX_12100_i707_i507
20.40363867
SLX_12100_i707_i508
<NA>
SLX_12100_i707_i517
1.537740298
SLX_12100_i708_i502
<NA>
SLX_12100_i708_i504
<NA>
SLX_12100_i708_i505
10.24219888
SLX_12100_i708_i506
61.65428491
SLX_12100_i708_i507
2.388567516
SLX_12100_i708_i508
1.11866837
SLX_12100_i708_i517
<NA>
SLX_12100_i709_i502
84.87846819
SLX_12100_i709_i503
49.39924658
SLX_12100_i709_i504
0.227605988
SLX_12100_i709_i505
37.82399987
SLX_12100_i709_i506
6.807672764
SLX_12100_i709_i507
0
SLX_12100_i709_i508
84.58015502
SLX_12100_i709_i517
<NA>
SLX_12100_i710_i502
76.92521644
SLX_12100_i710_i503
<NA>
SLX_12100_i710_i504
<NA>
SLX_12100_i710_i505
53.65480723
SLX_12100_i710_i506
22.70207271
SLX_12100_i710_i507
9.8468191
SLX_12100_i710_i508
51.33378868
SLX_12100_i710_i517
<NA>
SLX_12100_i711_i502
<NA>
SLX_12100_i711_i503
28.69816309
SLX_12100_i711_i504
<NA>
SLX_12100_i711_i506
1.656404451
SLX_12100_i711_i507
<NA>
SLX_12100_i711_i508
<NA>
SLX_12100_i711_i517
<NA>
SLX_12100_i712_i502
59.22134752
SLX_12100_i712_i503
52.81707104
SLX_12100_i712_i504
1.013338275
SLX_12100_i712_i505
1.000955539
SLX_12100_i712_i506
0.673295015
SLX_12100_i712_i507
13.90172358
SLX_12100_i712_i508
0.281432475
SLX_14251_i701_i502
25.50605057
SLX_14251_i701_i503
<NA>
SLX_14251_i701_i504
1.283897944
SLX_14251_i701_i505
<NA>
SLX_14251_i701_i506
9.034121294
SLX_14251_i701_i507
<NA>
SLX_14251_i701_i508
<NA>
SLX_14251_i701_i517
<NA>
SLX_14251_i702_i502
<NA>
SLX_14251_i702_i503
<NA>
SLX_14251_i702_i504
6.758001871
SLX_14251_i702_i505
1.109515561
SLX_14251_i702_i506
<NA>
SLX_14251_i702_i507
<NA>
SLX_14251_i702_i508
<NA>
SLX_14251_i702_i517
<NA>
SLX_14251_i703_i502
<NA>
SLX_14251_i703_i505
<NA>
SLX_14251_i703_i506
0.793321518
SLX_14251_i703_i507
<NA>
SLX_14251_i703_i508
29.0454459
SLX_14251_i703_i517
<NA>
SLX_14251_i704_i502
34.41450915
SLX_14251_i704_i503
1.088973394
SLX_14251_i704_i504
1.424888809
SLX_14251_i704_i505
0.825828051
SLX_14251_i704_i506
<NA>
SLX_14251_i704_i507
<NA>
SLX_14251_i704_i508
22.47583683
SLX_14251_i704_i517
19.99461322
SLX_14251_i705_i502
3.981433875
SLX_14251_i705_i503
1.382591289
SLX_14251_i705_i504
29.7646755
SLX_14251_i705_i505
22.50429502
SLX_14251_i705_i506
5.651177785
SLX_14251_i705_i507
<NA>
SLX_14251_i705_i508
23.25967007
SLX_14251_i705_i517
37.03714245
SLX_14251_i706_i502
1.590423549
SLX_14251_i706_i503
1.612251443
SLX_14251_i706_i504
0.845289381
SLX_14251_i706_i505
38.54100011
SLX_14251_i706_i506
0.568460427
SLX_14251_i706_i507
34.99047174
SLX_14251_i706_i508
0.797243913
SLX_14251_i706_i517
40.047026
SLX_14251_i707_i502
1.026614417
SLX_14251_i707_i503
3.454441829
SLX_14251_i707_i504
35.57231827
SLX_14251_i707_i505
21.11857844
SLX_14251_i707_i506
6.21855329
SLX_14251_i707_i507
1.351896085
SLX_14251_i707_i508
26.20585776
SLX_14251_i707_i517
39.31469375
SLX_14251_i708_i502
2.262331296
SLX_14251_i708_i503
27.59090282
SLX_14251_i708_i504
40.78698321
SLX_14251_i708_i505
35.29203109
SLX_14251_i708_i506
34.83781458
SLX_14251_i708_i507
40.10797005
SLX_14251_i708_i508
2.317748843
SLX_14251_i708_i517
40.09457648
SLX_14251_i709_i502
67.63240837
SLX_14251_i709_i503
86.0933035
SLX_14251_i709_i504
61.67196136
SLX_14251_i709_i505
85.53142528
SLX_14251_i709_i506
69.73795475
SLX_14251_i709_i507
66.07800116
SLX_14251_i709_i508
62.3983223
SLX_14251_i709_i517
86.18372624
SLX_14251_i710_i502
71.18272421
SLX_14251_i710_i503
58.48688407
SLX_14251_i710_i504
78.23350537
SLX_14251_i710_i505
85.08673304
SLX_14251_i710_i506
84.83061791
SLX_14251_i710_i507
65.42739471
SLX_14251_i710_i508
57.77809649
SLX_14251_i710_i517
59.36240419
SLX_14251_i711_i502
85.15301811
SLX_14251_i711_i503
86.15537326
SLX_14251_i711_i504
81.18650614
SLX_14251_i711_i505
76.85060805
SLX_14251_i711_i506
84.56941861
SLX_14251_i711_i507
66.07709325
SLX_14251_i711_i508
58.49121467
SLX_14251_i711_i517
76.22723714
SLX_14251_i712_i502
85.72440593
SLX_14251_i712_i503
69.0193979
SLX_14251_i712_i504
85.30660791
SLX_14251_i712_i505
53.60119161
SLX_14251_i712_i506
62.33815798
SLX_14251_i712_i507
85.91186274
SLX_14251_i712_i508
69.7431214
SLX_14251_i712_i517
59.28081613
SLX_14252_i701_i502
<NA>
SLX_14252_i701_i503
<NA>
SLX_14252_i701_i504
<NA>
SLX_14252_i701_i505
4.001223835
SLX_14252_i701_i506
0.3514962
SLX_14252_i701_i507
<NA>
SLX_14252_i701_i508
12.82396269
SLX_14252_i701_i517
<NA>
SLX_14252_i702_i502
<NA>
SLX_14252_i702_i503
<NA>
SLX_14252_i702_i504
0.378527623
SLX_14252_i702_i505
0.639130014
SLX_14252_i702_i506
<NA>
SLX_14252_i702_i507
<NA>
SLX_14252_i702_i508
...
SLX_14252_i702_i517
<NA>
SLX_14252_i703_i502
66.14144791
SLX_14252_i703_i503
<NA>
SLX_14252_i703_i504
<NA>
SLX_14252_i703_i505
<NA>
SLX_14252_i703_i506
<NA>
SLX_14252_i703_i507
<NA>
SLX_14252_i703_i508
<NA>
SLX_14252_i703_i517
1.185619541
SLX_14252_i704_i502
48.61774771
SLX_14252_i704_i503
72.65465605
SLX_14252_i704_i504
<NA>
SLX_14252_i704_i507
<NA>
SLX_14252_i704_i508
<NA>
SLX_14252_i704_i517
85.39882975
SLX_14252_i705_i502
<NA>
SLX_14252_i705_i503
<NA>
SLX_14252_i705_i504
<NA>
SLX_14252_i705_i505
<NA>
SLX_14252_i705_i506
<NA>
SLX_14252_i705_i507
<NA>
SLX_14252_i705_i508
<NA>
SLX_14252_i705_i517
<NA>
SLX_14252_i706_i502
1.068829038
SLX_14252_i706_i503
<NA>
SLX_14252_i706_i504
<NA>
SLX_14252_i706_i505
<NA>
SLX_14252_i706_i506
<NA>
SLX_14252_i706_i507
<NA>
SLX_14252_i706_i508
<NA>
SLX_14252_i706_i517
8.457272567
SLX_14252_i707_i502
<NA>
SLX_14252_i707_i503
<NA>
SLX_14252_i707_i504
<NA>
SLX_14252_i707_i505
4.494268855
SLX_14252_i707_i506
<NA>
SLX_14252_i707_i507
<NA>
SLX_14252_i707_i508
<NA>
SLX_14252_i707_i517
<NA>
SLX_14252_i708_i502
4.226597031
SLX_14252_i708_i503
<NA>
SLX_14252_i708_i504
<NA>
SLX_14252_i708_i505
<NA>
SLX_14252_i708_i506
86.04334115
SLX_14252_i708_i507
<NA>
SLX_14252_i708_i508
<NA>
SLX_14252_i708_i517
0.934120427
SLX_14252_i709_i502
<NA>
SLX_14252_i709_i503
21.08667017
SLX_14252_i709_i504
<NA>
SLX_14252_i709_i505
<NA>
SLX_14252_i709_i506
<NA>
SLX_14252_i709_i507
1.262182922
SLX_14252_i709_i508
<NA>
SLX_14252_i709_i517
<NA>
SLX_14252_i710_i502
<NA>
SLX_14252_i710_i503
<NA>
SLX_14252_i710_i504
<NA>
SLX_14252_i710_i505
<NA>
SLX_14252_i710_i506
<NA>
SLX_14252_i710_i507
<NA>
SLX_14252_i710_i508
<NA>
SLX_14252_i710_i517
<NA>
SLX_14252_i711_i502
1.638917613
SLX_14252_i711_i503
<NA>
SLX_14252_i711_i504
2.282926236
SLX_14252_i711_i505
<NA>
SLX_14252_i711_i506
<NA>
SLX_14252_i711_i507
57.49122462
SLX_14252_i711_i508
<NA>
SLX_14252_i711_i517
43.8667141
SLX_14252_i712_i502
<NA>
SLX_14252_i712_i503
<NA>
SLX_14252_i712_i504
<NA>
SLX_14252_i712_i505
<NA>
SLX_14252_i712_i506
<NA>
SLX_14252_i712_i507
<NA>
SLX_14252_i712_i508
<NA>
SLX_14252_i712_i517
6.700665965
SLX_14253_i701_i502
<NA>
SLX_14253_i701_i503
17.35305594
SLX_14253_i701_i504
<NA>
SLX_14253_i701_i505
<NA>
SLX_14253_i701_i506
<NA>
SLX_14253_i701_i507
<NA>
SLX_14253_i701_i508
<NA>
SLX_14253_i701_i517
<NA>
SLX_14253_i702_i502
52.82125401
SLX_14253_i702_i503
<NA>
SLX_14253_i702_i504
76.21090417
SLX_14253_i702_i505
4.520923862
SLX_14253_i702_i506
<NA>
SLX_14253_i702_i507
<NA>
SLX_14253_i702_i508
<NA>
SLX_14253_i702_i517
67.56291172
SLX_14253_i703_i502
<NA>
SLX_14253_i703_i503
<NA>
SLX_14253_i703_i504
<NA>
SLX_14253_i703_i505
<NA>
SLX_14253_i703_i506
<NA>
SLX_14253_i703_i507
85.70857599
SLX_14253_i703_i508
<NA>
SLX_14253_i703_i517
<NA>
SLX_14253_i704_i502
<NA>
SLX_14253_i704_i503
<NA>
SLX_14253_i704_i504
<NA>
SLX_14253_i704_i505
<NA>
SLX_14253_i704_i506
53.61734865
SLX_14253_i704_i507
<NA>
SLX_14253_i704_i508
<NA>
SLX_14253_i704_i517
<NA>
SLX_14253_i705_i502
<NA>
SLX_14253_i705_i503
<NA>
SLX_14253_i705_i504
<NA>
SLX_14253_i705_i505
<NA>
SLX_14253_i705_i506
<NA>
SLX_14253_i705_i507
<NA>
SLX_14253_i705_i508
<NA>
SLX_14253_i705_i517
1.715789222
SLX_14253_i706_i502
<NA>
SLX_14253_i706_i503
<NA>
SLX_14253_i706_i504
53.01526449
SLX_14253_i706_i505
<NA>
SLX_14253_i706_i506
<NA>
SLX_14253_i706_i507
<NA>
SLX_14253_i706_i508
31.17474744
SLX_14253_i706_i517
<NA>
SLX_14253_i707_i502
<NA>
SLX_14253_i707_i503
<NA>
SLX_14253_i707_i504
85.52909743
SLX_14253_i707_i505
85.84928458
SLX_14253_i707_i506
<NA>
SLX_14253_i707_i507
<NA>
SLX_14253_i707_i508
<NA>
SLX_14253_i707_i517
85.74228183
SLX_14253_i708_i502
49.68498932
SLX_14253_i708_i503
1.079153693
SLX_14253_i708_i504
<NA>
SLX_14253_i708_i505
<NA>
SLX_14253_i708_i506
<NA>
SLX_14253_i708_i507
<NA>
SLX_14253_i708_i508
20.01114307
SLX_14253_i708_i517
<NA>
SLX_14253_i709_i502
49.56628864
SLX_14253_i709_i503
<NA>
SLX_14253_i709_i504
<NA>
SLX_14253_i709_i505
13.51169843
SLX_14253_i709_i506
<NA>
SLX_14253_i709_i507
<NA>
SLX_14253_i709_i508
<NA>
SLX_14253_i709_i517
<NA>
SLX_14253_i710_i502
<NA>
SLX_14253_i710_i503
<NA>
SLX_14253_i710_i504
<NA>
SLX_14253_i710_i505
85.78794587
SLX_14253_i710_i506
<NA>
SLX_14253_i710_i507
<NA>
SLX_14253_i710_i508
<NA>
SLX_14253_i710_i517
<NA>
SLX_14253_i711_i502
<NA>
SLX_14253_i711_i503
52.82123089
SLX_14253_i711_i504
<NA>
SLX_14253_i711_i505
<NA>
SLX_14253_i711_i506
<NA>
SLX_14253_i711_i507
<NA>
SLX_14253_i711_i508
<NA>
SLX_14253_i711_i517
<NA>
SLX_14253_i712_i502
57.5232038
SLX_14253_i712_i503
<NA>
SLX_14253_i712_i504
<NA>
SLX_14253_i712_i505
<NA>
SLX_14253_i712_i506
<NA>
SLX_14253_i712_i507
<NA>
SLX_14253_i712_i508
<NA>
SLX_14253_i712_i517
<NA>
SLX_14254_i701_i502
<NA>
SLX_14254_i701_i503
<NA>
SLX_14254_i701_i504
1.763689235
SLX_14254_i701_i505
86.34070954
SLX_14254_i701_i506
<NA>
SLX_14254_i701_i507
<NA>
SLX_14254_i701_i508
57.51780473
SLX_14254_i701_i517
1.608888613
SLX_14254_i702_i502
<NA>
SLX_14254_i702_i503
<NA>
SLX_14254_i702_i504
<NA>
SLX_14254_i702_i505
<NA>
SLX_14254_i702_i506
64.65366714
SLX_14254_i702_i507
<NA>
SLX_14254_i702_i508
<NA>
SLX_14254_i702_i517
<NA>
SLX_14254_i703_i502
<NA>
SLX_14254_i703_i503
<NA>
SLX_14254_i703_i504
<NA>
SLX_14254_i703_i505
<NA>
SLX_14254_i703_i506
<NA>
SLX_14254_i703_i507
<NA>
SLX_14254_i703_i508
86.35105585
SLX_14254_i703_i517
<NA>
SLX_14254_i704_i502
4.253786819
In [47]:
seu$Pseudotime_PSE_rank <- rank(seu$Pseudotime_PSE)
In [48]:
FeaturePlot(seu, reduction = "umap", features = "Pseudotime_PSE_rank")
In [49]:
saveRDS(seu,"./tricycle/seu_Phloem_1375_cells.rds")
In [50]:
seu_PSE <- subset(seu, cells=colnames(seu)[!is.na(seu$Cluster_PSE)])
In [51]:
seu_PSE
An object of class Seurat 
32824 features across 758 samples within 1 assay 
Active assay: RNA (32824 features, 27218 variable features)
 2 dimensional reductions calculated: pca, umap
In [52]:
DimPlot(seu_PSE, reduction = "umap", group.by = "tricycleCCStage")
In [53]:
FeaturePlot(seu_PSE, reduction = "umap", features = "Pseudotime_PSE")
In [54]:
FeaturePlot(seu_PSE, reduction = "umap", features = "Pseudotime_PSE_rank")
In [55]:
saveRDS(seu_PSE,"./tricycle/seu_PSE_758_cells.rds")
In [56]:
## All cells (PSE, MSE, Procambium)
pltsg <- function(gene){
    keygene.idx <- which(rownames(seu@assays$RNA@data) == gene)
    fit.l <- fit_periodic_loess(seu$tricyclePosition*pi,
                            seu@assays$RNA@data[keygene.idx,],
                            plot = FALSE, span=0.3,
                       x_lab = "Cell cycle position \u03b8", y_lab = "log2(keygene)")
    fit <- fit.l$pred.df
    options(repr.plot.width=8, repr.plot.height=6)
    print(ggplot(fit, aes(x=x,y=y)) + geom_path(size = 1, alpha = 1)+ labs(x = "θ", y = "smoothed gene expression", title = gene) + 
                scale_x_continuous(limits = c(0, 2 * pi), breaks = c(0, 
                    pi/2, pi, 3 * pi/2, 2 * pi), labels = paste0(c(0, 
                    0.5, 1, 1.5, 2), "π"), name = "θ")+
      #              scale_x_continuous(limits = c(0, 200), breaks = c(0, 
      #              50, 100, 150, 200), labels = paste0(c(0, 
      #              0.5, 1, 1.5, 2), "π"), name = "θ")+
      theme_bw(base_size = 14)+
      theme(legend.title=element_blank()))
}
In [57]:
## DWF4
pltsg('AT3G50660')
In [58]:
## OPS
pltsg('AT3G09070')
In [59]:
## CPD
pltsg('AT5G05690')
In [60]:
## Only PSE
pltsg <- function(gene){
    keygene.idx <- which(rownames(seu_PSE@assays$RNA@data) == gene)
    fit.l <- fit_periodic_loess(seu_PSE$tricyclePosition*pi,
                            seu_PSE@assays$RNA@data[keygene.idx,],
                            plot = FALSE, span=0.3,
                       x_lab = "Cell cycle position \u03b8", y_lab = "log2(keygene)")
    fit <- fit.l$pred.df
    options(repr.plot.width=8, repr.plot.height=6)
    print(ggplot(fit, aes(x=x,y=y)) + geom_path(size = 1, alpha = 1)+ labs(x = "θ", y = "smoothed gene expression", title = gene) + 
                scale_x_continuous(limits = c(0, 2 * pi), breaks = c(0, 
                    pi/2, pi, 3 * pi/2, 2 * pi), labels = paste0(c(0, 
                    0.5, 1, 1.5, 2), "π"), name = "θ")+
      #              scale_x_continuous(limits = c(0, 200), breaks = c(0, 
      #              50, 100, 150, 200), labels = paste0(c(0, 
      #              0.5, 1, 1.5, 2), "π"), name = "θ")+
      theme_bw(base_size = 14)+
      theme(legend.title=element_blank()))
}
In [61]:
## DWF4
pltsg('AT3G50660')
In [62]:
## OPS
pltsg('AT3G09070')
In [63]:
## CPD
pltsg('AT5G05690')
In [64]:
seu_mer_PSE <- subset(seu_PSE, cells=colnames(seu_PSE)[which(seu_PSE$Pseudotime_PSE_rank <= 410)])
In [65]:
seu_mer_PSE
An object of class Seurat 
32824 features across 410 samples within 1 assay 
Active assay: RNA (32824 features, 27218 variable features)
 2 dimensional reductions calculated: pca, umap
In [66]:
saveRDS(seu_mer_PSE,"./tricycle/seu_mer_PSE_410_cells.rds")
In [67]:
## Only meristem PSE
pltsg <- function(gene){
    keygene.idx <- which(rownames(seu_mer_PSE@assays$RNA@data) == gene)
    fit.l <- fit_periodic_loess(seu_mer_PSE$tricyclePosition*pi,
                            seu_mer_PSE@assays$RNA@data[keygene.idx,],
                            plot = FALSE, span=0.3,
                       x_lab = "Cell cycle position \u03b8", y_lab = "log2(keygene)")
    fit <- fit.l$pred.df
    options(repr.plot.width=8, repr.plot.height=6)
    print(ggplot(fit, aes(x=x,y=y)) + geom_path(size = 1, alpha = 1)+ labs(x = "θ", y = "smoothed gene expression", title = gene) + 
                scale_x_continuous(limits = c(0, 2 * pi), breaks = c(0, 
                    pi/2, pi, 3 * pi/2, 2 * pi), labels = paste0(c(0, 
                    0.5, 1, 1.5, 2), "π"), name = "θ")+
      #              scale_x_continuous(limits = c(0, 200), breaks = c(0, 
      #              50, 100, 150, 200), labels = paste0(c(0, 
      #              0.5, 1, 1.5, 2), "π"), name = "θ")+
      theme_bw(base_size = 14)+
      theme(legend.title=element_blank()))
}
In [68]:
## DWF4
pltsg('AT3G50660')
In [69]:
## OPS
pltsg('AT3G09070')
In [70]:
## CPD
pltsg('AT5G05690')

Phloem Atlas (Nature Plants 2022, 10k cells)¶

In [71]:
sce <- readRDS("./scRNA-seq/Published_Data/phloem_atlas_sce.rds")
In [72]:
meta <- read.csv('./scRNA-seq/Published_Data/phloem_atlas_cell_metadata.csv')
In [73]:
head(meta)
A data.frame: 6 x 9
cell_idsamplebarcodetotal_countsdetected_genesclusterannotationUMAP1UMAP2
<chr><chr><chr><int><int><int><chr><dbl><dbl>
1APL_AAACCCAAGTAACGTA-1APLAAACCCAAGTAACGTA-13560052109early/undetermined 2.2210797-3.4197492
2APL_AAACCCACAACCAACT-1APLAAACCCACAACCAACT-1 661630582PSE 6.6928352 5.7139691
3APL_AAACCCACAGTTGAAA-1APLAAACCCACAGTTGAAA-11278136691MSE 0.9346761 0.6447128
4APL_AAACCCAGTGCCTTTC-1APLAAACCCAGTGCCTTTC-1 713831004PPP -5.9457399 1.5383996
5APL_AAACCCAGTTTCGACA-1APLAAACCCAGTTTCGACA-1 736425296PSE 5.8819741 8.5710071
6APL_AAACCCATCTGCGAGC-1APLAAACCCATCTGCGAGC-1 711428693CC -0.7623315-0.4224323
In [74]:
rownames(meta) <- meta$cell_id
In [75]:
nrow(meta)
10204
In [76]:
sce
class: SingleCellExperiment 
dim: 17396 10204 
metadata(4): nucgenes genevar hvgs hvgs_vst
assays(4): counts logcounts vst logvst
rownames(17396): AT1G01010 AT1G01020 ... AT5G67630 AT5G67640
rowData names(10): ID chrom ... go_mitosis go_phloem_dev
colnames(10204): APL_AAACCCAAGTAACGTA-1 APL_AAACCCACAACCAACT-1 ...
  sAPL_TTCGGTCAGCGACAGT-1 sAPL_TTTACCATCCTCTGCA-1
colData names(38): Sample Barcode ... slingPseudotime_4
  slingPseudotime_5
reducedDimNames(36): MNN_logcounts MNN_logvst ... DIFFMAP_MNN_logcounts
  DIFFMAP_MNN_logvst
mainExpName: NULL
altExpNames(0):
In [77]:
seu <- as.Seurat(sce)
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from MNN_logcounts_ to MNNlogcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to MNNlogcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from MNN_logvst_ to MNNlogvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to MNNlogvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from PC to PC_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to PC_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from PC to PC_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to PC_”
Warning message:
“Cannot add objects with duplicate keys (offending key: PC_), setting key to 'pca_logvst_'”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP7_logcounts_ to UMAP7logcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP7logcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP7_logvst_ to UMAP7logvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP7logvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP7_MNN_logcounts_ to UMAP7MNNlogcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP7MNNlogcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP7_MNN_logvst_ to UMAP7MNNlogvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP7MNNlogvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP15_logcounts_ to UMAP15logcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP15logcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP15_logvst_ to UMAP15logvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP15logvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP15_MNN_logcounts_ to UMAP15MNNlogcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP15MNNlogcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP15_MNN_logvst_ to UMAP15MNNlogvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP15MNNlogvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP30_logcounts_ to UMAP30logcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP30logcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP30_logvst_ to UMAP30logvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP30logvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP30_MNN_logcounts_ to UMAP30MNNlogcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP30MNNlogcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP30_MNN_logvst_ to UMAP30MNNlogvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP30MNNlogvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP100_logcounts_ to UMAP100logcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP100logcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP100_logvst_ to UMAP100logvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP100logvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP100_MNN_logcounts_ to UMAP100MNNlogcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP100MNNlogcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP100_MNN_logvst_ to UMAP100MNNlogvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP100MNNlogvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE15_logcounts_ to TSNE15logcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE15logcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE15_logvst_ to TSNE15logvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE15logvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE15_MNN_logcounts_ to TSNE15MNNlogcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE15MNNlogcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE15_MNN_logvst_ to TSNE15MNNlogvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE15MNNlogvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE30_logcounts_ to TSNE30logcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE30logcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE30_logvst_ to TSNE30logvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE30logvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE30_MNN_logcounts_ to TSNE30MNNlogcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE30MNNlogcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE30_MNN_logvst_ to TSNE30MNNlogvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE30MNNlogvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE60_logcounts_ to TSNE60logcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE60logcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE60_logvst_ to TSNE60logvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE60logvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE60_MNN_logcounts_ to TSNE60MNNlogcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE60MNNlogcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE60_MNN_logvst_ to TSNE60MNNlogvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE60MNNlogvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from DC to DC_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to DC_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from DC to DC_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to DC_”
Warning message:
“Cannot add objects with duplicate keys (offending key: DC_), setting key to 'diffmap_logvst_'”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from DC to DC_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to DC_”
Warning message:
“Cannot add objects with duplicate keys (offending key: DC_), setting key to 'diffmap_mnn_logcounts_'”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from DC to DC_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to DC_”
Warning message:
“Cannot add objects with duplicate keys (offending key: DC_), setting key to 'diffmap_mnn_logvst_'”
In [78]:
seu
An object of class Seurat 
17396 features across 10204 samples within 1 assay 
Active assay: originalexp (17396 features, 0 variable features)
 36 dimensional reductions calculated: MNN_logcounts, MNN_logvst, PCA_logcounts, PCA_logvst, UMAP7_logcounts, UMAP7_logvst, UMAP7_MNN_logcounts, UMAP7_MNN_logvst, UMAP15_logcounts, UMAP15_logvst, UMAP15_MNN_logcounts, UMAP15_MNN_logvst, UMAP30_logcounts, UMAP30_logvst, UMAP30_MNN_logcounts, UMAP30_MNN_logvst, UMAP100_logcounts, UMAP100_logvst, UMAP100_MNN_logcounts, UMAP100_MNN_logvst, TSNE15_logcounts, TSNE15_logvst, TSNE15_MNN_logcounts, TSNE15_MNN_logvst, TSNE30_logcounts, TSNE30_logvst, TSNE30_MNN_logcounts, TSNE30_MNN_logvst, TSNE60_logcounts, TSNE60_logvst, TSNE60_MNN_logcounts, TSNE60_MNN_logvst, DIFFMAP_logcounts, DIFFMAP_logvst, DIFFMAP_MNN_logcounts, DIFFMAP_MNN_logvst
In [79]:
saveRDS(seu,"./tricycle/phloem_atlas_seu4.rds")
In [80]:
seu <- CreateSeuratObject(seu@assays$originalexp@counts, meta.data=meta)
In [81]:
seu <- NormalizeData(seu)
seu <- FindVariableFeatures(seu, selection.method = "vst", nfeatures = nrow(seu))
all.genes <- rownames(seu)
seu <- ScaleData(seu, features = all.genes)
seu <- RunPCA(seu,features = VariableFeatures(object = seu))
Centering and scaling data matrix

PC_ 1 
Positive:  AT5G16130, AT1G27400, AT2G37190, AT3G45030, AT5G27850, AT4G29410, AT2G34480, AT3G16080, AT3G62870, AT5G22440 
	   AT1G22780, AT2G44120, AT3G11940, AT5G15200, AT1G09690, AT3G61110, AT5G20290, AT4G36130, AT2G39460, AT2G01250 
	   AT3G22230, AT5G59850, AT3G07110, AT5G60670, AT5G48760, AT1G69620, AT5G52650, AT5G39740, AT1G09590, AT5G45775 
Negative:  AT1G76180, AT1G62480, AT5G65020, AT1G65930, AT2G02130, AT1G49240, AT5G56540, AT2G19760, AT3G21770, AT3G13520 
	   AT1G10682, AT2G45820, AT1G20440, AT3G14310, AT1G05850, AT3G53420, AT5G17920, AT3G05890, AT1G12080, AT2G44790 
	   AT2G30870, AT5G11970, AT3G23820, AT3G11930, AT3G56240, AT5G28050, AT1G35720, AT2G47730, AT1G55330, AT4G26320 
PC_ 2 
Positive:  AT1G29470, AT3G58120, AT5G23860, AT2G45070, AT1G23490, AT5G14030, AT2G27030, AT1G70490, AT2G47170, AT5G15350 
	   AT3G12610, AT5G48580, AT1G78040, AT1G04680, AT1G04820, AT3G03160, AT1G10630, AT2G27130, AT2G21160, AT3G43810 
	   AT1G50010, AT5G17190, AT4G14960, AT4G24920, AT3G54810, AT2G19460, AT4G32460, AT5G03345, AT4G24190, AT3G25220 
Negative:  AT4G14010, AT2G37040, AT2G40890, AT4G34600, AT3G21240, AT5G57685, AT1G32100, AT5G64260, AT2G44300, AT5G54160 
	   AT4G34050, AT5G48930, AT3G27890, AT1G44800, AT2G36830, AT2G30860, AT1G27030, AT2G19110, AT2G40900, AT3G23090 
	   AT1G73660, AT1G18140, AT2G30490, AT3G26520, AT1G48320, AT3G10340, AT1G16350, AT1G06550, AT4G15340, AT1G48750 
PC_ 3 
Positive:  AT3G28860, AT5G19780, AT4G19410, AT1G15690, AT1G09640, AT5G60390, AT5G02500, AT3G19820, AT1G07940, AT2G31360 
	   AT4G13940, AT2G36530, AT4G38680, AT2G14890, AT4G34200, AT5G08690, AT5G14040, AT5G44340, AT3G20050, AT3G11830 
	   AT5G62890, AT5G08670, AT3G04120, AT2G34250, AT1G04430, AT1G20330, AT3G13460, AT1G24510, AT5G20890, AT5G46290 
Negative:  AT4G24960, AT1G62045, AT3G17730, AT5G18130, AT2G35585, AT5G14110, AT1G69970, AT3G49360, AT5G50290, AT3G54040 
	   AT5G02140, AT3G01940, AT1G33475, AT2G18380, AT1G18210, AT2G27740, AT1G06490, AT4G17510, AT1G65910, AT3G02677 
	   AT3G18670, AT1G63310, AT1G11570, AT1G75720, AT3G08690, AT4G05220, AT1G05610, AT1G31935, AT3G03170, AT5G01070 
PC_ 4 
Positive:  AT1G22840, AT4G32470, AT1G31812, AT1G72020, AT1G51650, AT4G21105, AT2G31490, AT5G64350, AT4G29480, AT1G15120 
	   AT3G46430, AT5G47030, AT5G53560, AT4G20150, AT1G27450, AT3G52730, AT4G16450, AT5G08060, AT3G52300, AT4G34700 
	   AT5G47890, AT4G37830, AT1G05205, AT5G59613, AT4G30010, AT1G76200, AT1G14450, AT3G10860, AT4G09030, AT5G05370 
Negative:  AT3G55760, AT3G24240, AT1G10760, AT5G14640, AT4G30020, AT2G37590, AT1G31880, AT1G06490, AT5G19090, AT1G65910 
	   AT1G55760, AT1G54030, AT2G17840, AT5G61960, AT3G02260, AT4G05430, AT3G54040, AT5G01890, AT3G49360, AT5G60200 
	   AT1G63300, AT2G18380, AT1G64390, AT4G39400, AT1G05200, AT1G05610, AT3G41768, AT1G72700, AT5G14110, AT2G37080 
PC_ 5 
Positive:  AT5G26260, AT3G16420, AT1G54000, AT2G41800, AT2G47140, AT3G53480, AT2G43610, AT3G13790, AT5G60950, AT5G10130 
	   AT3G15950, AT3G09260, AT5G02260, AT3G20370, AT1G53180, AT3G16460, AT2G04160, AT5G26280, AT2G17280, AT2G43535 
	   AT2G30340, AT1G66270, AT1G78340, AT5G35940, AT1G54010, AT5G01040, AT5G14730, AT3G29034, AT4G23590, AT2G45400 
Negative:  AT2G26730, AT4G20270, AT1G68430, AT1G60390, AT3G15680, AT5G60200, AT2G16850, AT4G22120, AT1G72150, AT2G45470 
	   AT3G52490, AT1G80690, AT5G25490, AT3G04670, AT5G57110, AT1G73850, AT3G63300, AT1G05470, AT2G18380, AT1G65910 
	   AT3G28455, AT2G01910, AT5G54660, AT4G36620, AT1G70670, AT1G56230, AT1G63300, AT2G44830, AT3G45610, AT3G59680 

In [82]:
seu <- FindNeighbors(seu, dims = 1:10)
seu <- FindClusters(seu, resolution = 0.5)
Computing nearest neighbor graph

Computing SNN

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 10204
Number of edges: 303357

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9067
Number of communities: 14
Elapsed time: 1 seconds
In [83]:
seu <- RunUMAP(seu, dims = 1:10)
10:59:20 UMAP embedding parameters a = 0.9922 b = 1.112

10:59:21 Read 10204 rows and found 10 numeric columns

10:59:21 Using Annoy for neighbor search, n_neighbors = 30

10:59:21 Building Annoy index with metric = cosine, n_trees = 50

0%   10   20   30   40   50   60   70   80   90   100%

[----|----|----|----|----|----|----|----|----|----|

*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
|

10:59:22 Writing NN index file to temp file /tmp/RtmpSzDebc/file16ca8b2cbb33de

10:59:22 Searching Annoy index using 1 thread, search_k = 3000

10:59:25 Annoy recall = 100%

10:59:27 Commencing smooth kNN distance calibration using 1 thread
 with target n_neighbors = 30

10:59:30 Initializing from normalized Laplacian + noise (using RSpectra)

10:59:30 Commencing optimization for 200 epochs, with 388778 positive edges

10:59:37 Optimization finished

In [84]:
head(seu@meta.data)
A data.frame: 6 x 14
orig.identnCount_RNAnFeature_RNAcell_idsamplebarcodetotal_countsdetected_genesclusterannotationUMAP1UMAP2RNA_snn_res.0.5seurat_clusters
<fct><dbl><int><chr><chr><chr><int><int><int><chr><dbl><dbl><fct><fct>
APL_AAACCCAAGTAACGTA-1APL353835176APL_AAACCCAAGTAACGTA-1APLAAACCCAAGTAACGTA-13560052109early/undetermined 2.2210797-3.41974928 8
APL_AAACCCACAACCAACT-1APL 64563020APL_AAACCCACAACCAACT-1APLAAACCCACAACCAACT-1 661630582PSE 6.6928352 5.71396919 9
APL_AAACCCACAGTTGAAA-1APL127313650APL_AAACCCACAGTTGAAA-1APLAAACCCACAGTTGAAA-11278136691MSE 0.9346761 0.64471288 8
APL_AAACCCAGTGCCTTTC-1APL 69543072APL_AAACCCAGTGCCTTTC-1APLAAACCCAGTGCCTTTC-1 713831004PPP -5.9457399 1.53839964 4
APL_AAACCCAGTTTCGACA-1APL 72762483APL_AAACCCAGTTTCGACA-1APLAAACCCAGTTTCGACA-1 736425296PSE 5.8819741 8.57100711111
APL_AAACCCATCTGCGAGC-1APL 70452849APL_AAACCCATCTGCGAGC-1APLAAACCCATCTGCGAGC-1 711428693CC -0.7623315-0.42243231 1
In [85]:
DimPlot(seu, reduction = "umap", group.by = "annotation")
In [86]:
sce <- as.SingleCellExperiment(seu)

## Tricycle estimation
sce <- project_cycle_space(sce,ref.m = new.ref)
sce <- estimate_cycle_position(sce)
The number of projection genes found in the new data is 242.

In [87]:
sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
    sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
    sce$tricyclePosition <- sce$tricyclePosition*pi
    
    ## Plot
    print(plot_ccposition_den(sce$tricyclePosition,
                    sce$orig.ident, 'sample',
                    bw = 10, fig.title = "Kernel density of \u03b8") +
      theme_bw(base_size = 14))
    
    ## Assign position to stage
    tricycleCCStage <- sce$tricyclePosition/pi

    G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
    S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
    G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
                    
    tricycleCCStage[G1_idx]="G1/G0"
    tricycleCCStage[S_idx]="S"
    tricycleCCStage[G2_idx]="G2M"
    
    sce$tricycleCCStage <- tricycleCCStage
    
    ## Schwabe Stage
    #sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
    
    ## Table
    #seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
    #print(table(seu$time.anno.Li, sce$tricycleCCStage))
    
    ## Prepare seu
    #seu <- prep(seu)
    #plot_anno(seu)
    
    ## Sce to seu
    seu$tricyclePosition <- sce$tricyclePosition/pi
    seu$tricycleCCStage <- sce$tricycleCCStage
    #seu$SchwabeCCStage <- sce$CCStage
In [88]:
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                       point.size = 4, point.alpha = 1) +
  theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 1)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
In [89]:
DimPlot(seu, reduction = "umap", group.by = "tricycleCCStage")
In [90]:
saveRDS(seu,"./tricycle/phloem_atlas_seu4_simplified.rds")

DE analysis¶

In [ ]:
# Set the identity of interest to cell type
Idents(rc.integrated) <- "tricycleCCStage"
# Use SCT assay 
DefaultAssay(rc.integrated) <- "integrated"
de.list <- FindAllMarkers(rc.integrated,
                            max.cells.per.ident = 10000,
                            only.pos=T, 
                           test.use="wilcox")
Calculating cluster G2M

In [ ]:
de.list <- de.list %>% filter(p_val_adj < 0.01 & avg_log2FC > 0.25)
In [18]:
table(de.list$cluster)
  G2M G1/G0     S 
  133     9   103 
In [19]:
de.list %>% filter(cluster=="G1/G0")
A data.frame: 9 x 7
p_valavg_log2FCpct.1pct.2p_val_adjclustergene
<dbl><dbl><dbl><dbl><dbl><fct><chr>
AT5G028201.728994e-1202.70873320.1380.3784.236034e-118G1/G0AT5G02820
AT1G04710 2.012929e-902.83750440.1200.281 4.931675e-88G1/G0AT1G04710
AT2G03670 4.696481e-504.88834490.1380.325 1.150638e-47G1/G0AT2G03670
AT3G13060 2.300513e-251.98374730.1330.147 5.636256e-23G1/G0AT3G13060
AT2G42190 1.941702e-210.85646820.1260.202 4.757170e-19G1/G0AT2G42190
AT4G14980 3.124803e-082.96060000.1110.129 7.655767e-06G1/G0AT4G14980
AT1G05520 1.283976e-071.87092910.1830.239 3.145741e-05G1/G0AT1G05520
AT1G128001 1.191791e-060.72927460.1410.264 2.919888e-04G1/G0AT1G12800
AT5G63120 2.330971e-064.86042930.1100.179 5.710880e-04G1/G0AT5G63120
In [ ]: